/M2- ZVW 


THE USE OF LINEAR PROGRAMMING TECHNIQUES 
TO DESIGN OPTIMAL DIGITAL FILTERS FOR 
PULSE SHAPING AND CHANNEL EQUALIZATION 


. ’iC' 

CO 

fev. 

15 

t' i 

EX 

iOUTS 



Project Director- 


DONALD W. BURLAGE 
Research Associate 


April, 1972 


TECHNICAL REPORT NUMBER 142-102 


COMMUNICATION SYSTEMS GROUP 


BUREAU OF ENGINEERING RESEARCH 
UNIVERSITY OF ALABAMA UNIVERSITY, ALABAMA 


THE USE OF LINEAR PROGRAMMING TECHNIQUES 


TO DESIGN OPTIMAL DIGITAL FILTERS FOR 
PULSE SHAPING AND CHANNEL EQUALIZATION 


by 

RONALD C. HOUTS 
Project Director 

and 

DONALD W. BURLAGE 
Research Associate 


APRIL, 1972 

TECHNICAL REPORT NUMBER 142-102 


Prepared for 

National Aeronautics and Space Administration 
Marshall Space Flight Center 
Huntsville, Alabama 


Under Contract Number 
NAS8-20172 


Communication Systems Group 
Bureau of Engineering Research 
University of Alabama 



ABSTRACT 


A time-domain technique is developed to design finite-duration 
impulse-response digital filters using linear programming. Two 
related applications of this technique in data transmission systems 
are considered. The first is the design of pulse-shaping digital 
filters to generate or detect signaling waveforms transmitted over 
bandlimited channels that are assumed to have ideal low-pass or band- 
pass characteristics. The second is the design of digital filters to 
be used as preset equalizers in cascade with channels that have known 
impulse-response characteristics . 

A straightforward design procedure which can be used for both 
applications is established. Specifications for the design are 
expressed in terms of the desired impulse response output from either 
a pulse-shaping filter or an equalized transmission channel. Equa- 
tions which are linear functions of the unknown filter multiplier coef- 
ficients are then derived corresponding to the appropriate impulse 
response. Finally, these equations are constrained to have the speci- 
fied values at appropriate critical time points so that a linear pro- 
gramming routine can be applied to solve for the multiplier coeffi- 
cients . 

Example designs are presented which illustrate that excellent 
waveforms can be generated with frequency-sampling filters and the 
ease with which digital transversal filters can be designed for preset 
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equalization. Other results are included to depict the relation- 
ships and tradeoffs between such digital filter parameters as sampling 
interval and number of multipliers required to meet such equalization 
specifications as amount of residual intersymbol interference which 
can be tolerated and input energy which must be transmitted. 
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CHAPTER 1 


INTRODUCTION 

Recent developments in digital-circuit technology have indicated 
that the processing of signals with digital hardware, i.e., digital 
filtering, has some very significant advantages over the conventional 
analog signal processor; particularly, greater flexibility and accu- 
racy, smaller size, and potentially lower implementation cost [1], 
These advantages have resulted in a great upsurge in the research 
and development of digital systems to perform functions previously 
implemented only with analog hardware, and even more important, 
functions that were considered impossible due to the limitations of 
analog components. Emphasis has been placed on the design of soft- 
ware as well as digital hardware. Algorithms are being developed 
and programmed on digital computers to perform operations that pre- 
viously required complete systems of conventional analog design. 

This study will concentrate on one of these areas of research 
that is becoming increasingly important, namely the application of 
digital filters to solve problems of baseband data communications. 
Transmission of pulse-amplitude-modulated (PAM) signals over any 
realistic channel is impeded both by additive noise and time dis- 
persion resulting from the channel's nonlinear phase characteristic. 
Consequently, emphasis will be given to the design and implementation 
of pulse-shaping digital filters for the purpose of reducing these 
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two effects. Although filters have been designed to compensate for 
these two effects, most designs have required a complicated analog 
implementation. Hence, the intent of this study is to develop a 
simple-to-use algorithm for the time-domain design. of digital filters 
so that the resulting filter can be implemented easily and will have 
the desired compensating effect on a data transmission channel. 

A background for the study is provided in Section 1.1 by the 
discussion of problems inherent in baseband data transmission and 
the review in Section 1.2 of previous efforts to solve these prob- 
lems. The content of this study is outlined in Section 1.3. 

1.1 BASEBAND TRANSMISSION PROBLEMS 

The transmitted signal for PAM waveforms is given by the expres- 
sion 

s(t) = a_. k g(t - kt b ) , (1.1) 

k 

where the coefficient a^ k assumes one of M possible values, denoted 

th 

by the subscript j = 1,2,...,M, for the k interval. The time inter- 
val between transmitted pulses is t^, and the set of coefficients 

(a., }, which are called "input symbols", are generated by a trans- 
it 

mitter data source at the rate of 1/ t^ symbols per second. For 

example, transmission of binary data would require a coefficient set 

of at least two values, e.g., a,. = A, a„, = -A. 

Ik 2k 

This investigation will emphasize the use of digital filters 
for the generation and detection of the individual pulses g(t) that 
make up the transmitted signal. For example, if there is some sort 
of encoder-decoder device in the communication system, so as to 
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increase transmission efficiency, the symbols {a., } will be defined 

Jk 

to be the encoder output. The pulses will be assumed to have a base- 
band frequency spectrum, i.e., frequency components are present near, 
but not necessarily including, zero frequency. Any frequency trans- 
lation operation such as amplitude or frequency modulation of a car- 
rier will be assumed to be part of the transmission channel. 


1.1.1 Intersymbol Interference 

A fundamental impairment to the transmission of PAM waveforms 
is the frequency distortion and bandwidth limitations of realistic 
channels, e.g. , consider the transmitted and received pulses sketched 
in Fig. l-l(a). For this system, g(t) is a rectangular pulse of 
unit magnitude and width t^ seconds. The transmitted signal would 
thus consist of a sum of these pulses, each multiplied by a^. How- 
ever, as indicated, a channel with impulse response c(t) produces a 
time . spread for each pulse due to phase shifting of the spectral 
components. Each received pulse y(t) results from the convolution 
of g(t) and c(t), denoted hereafter as g(t) * c(t). Therefore, the 
received waveform r(t) consists of overlapping pulses, i.e., 


r(t) = a^ k y (t - kt fe ) . (1.2) 

k 

This time overlap results in a distortion called intersymbol inter- 
ference (ISI) which can cause errors to be made in the detection of 
pulses, even in the absence of noise in the channel. 

If the input symbols are to be determined by sampling the 
received signal at intervals of t^ seconds and comparing the resulting 
samples to threshold values, as is often done in simple PAM receivers. 




(c) Block Diagram of Noise-Free Baseband Systea 


Fig. 1-1. Baseband Data Transmission. 
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the nonzero magnitudes at ± t^ and 2t^ in Fig. 1-1 (a) represent ISI 
distortion introduced in the sample values of adjacent pulses. In 
this manner, the total peak ISI distortion caused by the pulse, 
which appears distributed in all other sample values, is defined as 

CO 

D k = Z a jk |yat b - kt b } i • (i - 3) 

£=-“ 

£ 4 k 

The same pulse waveshape is used to transmit all symbols over the 
same channel, so the distortion value 

\ 

D = f- (1.4) 

a jk 

is independent of the symbol being transmitted. Hence, reduction of 
intersymbol interference is equivalent to reducing the value of D. 

Complete elimination of ISI distortion, i.e., D = 0, would 
require a waveform y(t) at the channel output with zero magnitude at 
the sampling times of all other pulses, i.e., at integer multiples of 
t^ seconds measured with reference to the pulse peak. Such a situa- 
tion is illustrated by the overlapping pulses of Fig. l-l(b). These 
pulses, which correspond to the sequence of transmitted symbols -A, 

A, and A, must be summed to give the actual signal r(t) at the channel 
output. However, it is obvious from Fig. 1-1 (b) that the magnitude of 
r(t) would be identically -A, A, and A at the respective receiver 
sampling instants -t^, 0, and t^. Thus, the receiver decisions of the 
transmitted symbols, where the detector output is denoted by a^, are 
error free in the absence of channel noise. The received pulse shape 
of Fig. l-l(b) is called a raised-cosine pulse and is of a class that 
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Nyquist [2] showed would have zero intersymbol interference and require 

a bandwidth of approximately l/(2t, ) Hz. Pulses of this class are used 

b 

as optimum design goals and can only be approximated in practice. Fur- 
thermore, it is difficult to approach very closely this type of pulse 
with analog filters; however, excellent results can be obtained with 
digital filters. 

A technique will be developed to design pulse-shaping digital- 
filters that can be used to reduce the distortion D, As shown in Fig. 
1-1 (c) , the pulse g(t) is generated with a digital transmitter filter 
and applied to a known channel so that the channel output y(t) will 
closely approximate the raised-cosine pulse. Note that the k subscript 
for a . k and a^ as used in Fig. 1-1 (b) is dropped when the transmission 
of a single pulse is being considered. The analog equivalent of. the 
transmitter filter has been called an "equalizer" rather than a pulse- 
shaping filter, since it is usually thought of as a device which alters 
the composite frequency characteristic in order to reduce I§I distor- 
tion. For some applications the channel characteristic c(t) will be 
unknown, in which case it will be assumed that the channel is distor- 
tion free. Consequently, the pulse applied to the channel, rather than 
the channel output, is specified to be of the Nyquist class. A tech- 
nique for the design of pulse-shaping digital filters for this applica- 
tion will be developed in Chapter 3 and the design technique using 
known channel specifications in Chapter 4. 


1.1.2 Channel Noise 

The received signal as defined in (1.2) was assumed to be noise 
free. For many channels which employ coaxial cable or shielded, 
twisted pairs of wires, the signal-to-noise ratio is high and the 
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noise-free assumption is valid. Under these circumstances reduction 
of intersymbol interference is of prime importance. However, for low 
signal-to-noise ratios the emphasis must be placed on optimum detection 
of the transmitted pulse in the presence of noise. The channel noise 
is usually modeled by adding a term n(t) to the received signal, i.e., 

r (t) = ^ a k y (t - kt b > + n(t) , (1.5) 

k 

where n(t) is assumed to be a white Gaussian noise signal. 

It has been shown [3] that a matched-filter detector will give 
the optimum (minimum probability of detector error) reception of sig- 
nals in an additive white Gaussian noise background. The block diagram 
of such a detector is given in Fig. 1-2 for the detection of two 
equally likely equal-energy signals, y^(t) andy 2 (t). This detector 
includes two filters, each of which is "matched" to one of the possible 
signals, i.e., the filter unit-impulse responses are 

hi (t) = Y i (t b - t) i = 1,2 . (1.6) 

A sampler and decision circuit form the remaining portions of the detec- 
tor. During each time interval of t b seconds r(t) will consist of the 
additive noise plus either y^(t) or y 2 (t). At the end of each interval, 
the filter outputs are sampled to obtain the decision statistics and. 
S 2 . The best estimate of which signal is present, i.e., which symbol 
has been transmitted, is then obtained by comparing the two statistics 
as indicated in the block diagram. 

The matched-filter detector can be simplified further for signal 
sets where all possible signals are multiples of the same unique 
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function. For example, if it is assumed that the two signals imbedded 
in the noise are y^(t) = A y(t) and y^(t) = -A y(t), where the function 
y(t) is assumed to be identically zero outside the interval 0 <_ t <_ t^, 

then S, = — S_ at time t = t, . Hence, the decision rule can be modi- 

l 2 b 

fied so that all that is required is to determine whether is posi- 
tive or negative, which eliminates the need for one of the matched 
filters. The PAM baseband received signal (1.5) also consists of 
multiples of a single function y(t) plus the noise interference. 

Hence, if y(t) satisfies the time-limited constraint, the optimum 
receiver for detection of the PAM pulses is a filter, sampler, and 
threshold detector as illustrated in Fig. 1-3. Of course, it is 
impossible for the function y(t) to be truly time limited at the out- 
put of a bandlimited channel. Nevertheless, optimum detection of PAM 
baseband signals can be approximated with this type of matched-filter 
technique. The same technique used to design digital filters for 
transmission of pulses over unspecified channels will also be used to 
design digital filters to implement the matched-filter portion of 
this receiver. 

1.2 HISTORICAL REVIEW 

Efforts to improve the quality of PAM data transmission over 
the baseband channel range from Nyquist's work [4] on telegraph trans- 
mission theory in the 1920 's to the present day efforts to design all 
digital PAM transmitters and receivers. Summaries of these studies 
have been published [5,6] and some of the more pertinent results, 
namely theoretical developments in PAM data transmission, digitally 
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synthesized PAM systems, coding-decoding techniques, and adaptive 
equalization will be reviewed in this section. 

1.2.1 Theoretical Developments 

In 1928 Nyquist published a classic paper [2] which has formed 
the basis for the study of PAM data transmission. He developed crite- 
ria for the required overall frequency characteristic of a finite- 
bandwidth transmission channel so that distortionless transmission 
could be achieved. He also illustrated that the minimum bandwidth 
required for distortionless baseband data transmission is approxi- 
mately equal to the reciprocal of twice the symbol interval, i.e., 
l/2t^. The time interval, t^, is called the "Nyquist interval" and 
the corresponding bandwidth, l/2t^, the "Nyquist bandwidth". However, 
Nyquist' s technique required the solution of complicated Fourier inte- 
grals; consequently, an approximation technique called the "method of 
paired echoes" was developed by Wheeler [7]. He showed that small 
sinusoidal perturbations in either the channel amplitude or phase 
characteristic would result in echoes in the channel impulse-response; 
consequently, the location and amplitude of the echoes could be pre- 
dicted by inspection of either frequency characteristic. In 1954 
Sunde published a comprehensive summary [8,9] of the theoretical work 
of these early researchers. 

More recently, emphasis has been placed on techniques to reduce 
intersymbol interference by placing an equalizing device at the trans- 
mitter and/or receiver, or even in the center of the transmission chan 
nel. Bogert [10] developed what he called a "time reversal technique" 
which in theory would eliminate the phase distortion effects of any 
bisectional channel. He illustrated that the channel output could be 
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made independent of channel phase characteristics by placing an 
equalizer at the channel midpoint to reverse the time base of signal 
segments. However, the impractical technique he used for time rever- 
sal was to record signal segments at the channel midpoint on magnetic 
tape, and then run the recorded signal backwards past the reproduce 
heads during playback into the second half of the channel. 

Gerst & Diamond [11] developed a technique to specify the 
required impulse response of a transmitter filter so that the out- 
put of a known channel would be time-limited to the Nyquist interval 
of t^ seconds. Consequently, ISI distortion would be completely 
eliminated. Unfortunately, their technique is only applicable to 
channels that are not bandlimited, and requires maximum energy trans- 
fer at frequencies where the channel attenuation is greatest. 

Aein & Hancock [12] considered the problems of both background 
channel noise and ISI distortion. They assumed that the received 
pulse would only overlap in the following Nyquist interval, and then 
derived the impulse response for the matched-filter receiver which 
would be optimum under this assumption. Their calculated probability 
of-error curves illustrated the performance improvement of their 
receiver over the standard matched-filter detector of Fig. 1-3. 

George [13] removed the restriction on the location of ISI distor- 
tion and performed a similar analysis to define an optimum matched- 
filter receiver. He assumed the noise to be a sum of ISI distortion 
and channel noise. 

Gibby & Smith [14] extended the distortionless transmission 
criteria of Nyquist by removing the bandwidth restrictions. Their 
results are generalized to include channels with either no low 
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frequency response or bandwidths that exceed the Nyquist bandwidth. 
Constraints were derived on the real and imaginary parts of the chan- 
nel transfer characteristic to ensure distortionless transmission. 

The results of a number of efforts to jointly optimize both the 
transmitter and receiver filters have been published. Tufts [15] 
studied the problem of joint transmitter and receiver filter optimiza 
tion to minimize the mean-square-difference between the channel input 
symbols and the receiver output estimates for a finite sequence of 
data. The difference was assumed to result from the combined effects 
of intersymbol interference and channel noise. He showed that the 
optimum receiver was the cascade of a matched filter and a tapped 
delay line with time-varying coefficients. The corresponding opti- 
mum transmitted pulse waveform was required to be time limited to 
one Nyquist interval and to be a time-varying function. Aaron & 

Tufts [16] proved that this structure for minimizing the mean-square- 
difference also minimizes the average probability of detector error 
over a finite sequence of symbols. Coefficients of the tapped delay 
line calculated by both minimization techniques agree closely. How- 
ever, calculation of these tap coefficients requires the difficult 
solution of simultaneous nonlinear equations, and thus severely 
limits the technique for practical system design. Smith [17] removed 
the restriction that the transmitted pulse be time limited and also 
obtained a transmitter-receiver filter combination which was time 
invariant. Furthermore, his approach did not require the solution 
of nonlinear equations, but he could only minimize the mean-square- 
difference to an approximation. Berger & Tufts [18] extended the 
optimization study to include the effects of timing jitter and also 
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compared performance characteristics of the jointly optimized PAM 
systems to the theoretically attainable channel capacity as derived 
by Shannon. 

1.2.2 Digital Synthesis Developments 

The results of theoretical attempts to optimize PAM system 
performance as described in the previous subsection usually require 
the accurate generation of complicated signaling waveforms and the 
equally accurate synthesis of matched filters with specified impulse 
responses. Technical and economical problems associated with con- 
ventional analog techniques for waveform generation and matched filter 
realization have resulted in the recent attempts to design digital PAM 
transmitters and receivers. 

Voelcker [19] was one of the first to illustrate the use of 
digital hardware to generate PAM signals with his binary transversal 
filter (BTF) . His technique, which was actually a hybrid realization, 
was to replace the tapped delay line of the conventional transversal 
filter with a shift register and to use resistive summing networks to 
form the tap coefficients. He didn't consider the design problem of 
specifying filter coefficients but did analyze in detail the sources 
of distortion inherent in the BTF and how to reduce their effects. 
Since the summing portion of the filter is analog, the problems of 
resistor precision, loading effects, and voltage variations in the 
shift-register stages are still disadvantages to this type of signal 
generation. 

Van Gerwen & Van Der Wurf [20] incorporated the BTF and a digital 
modulator into an experimental vestigial-sideband data transmitter. 
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By using digital circuitry, they were able to implement the trans- 
mitter entirely with resistors and transistors. Consequently, the 

complete transmitter, consisting of 203 transistors and 172 resistors, 

2 

was integrated into one 2.1 x 2.7 mm chip. Satisfactory transmission 
was obtained up to a signaling rate of 1 Mbps. They also discussed 
the digital implementation of other modulation techniques such as 
frequency-shift keying and orthogonal modulation as well as the use 
of digital filters in receivers. 

Nowak & Schmid [21] studied the problem of designing nonrecursive 
digital filters* to have the raised-cosine frequency characteristic 
specified by Nyquist. They used a Fourier series expansion of the 
raised-cosine magnitude characteristic to determine the nonrecursive 
filter coefficients. A major portion of their study was devoted to 
measuring the filter's minimum stop-band attenuation and impulse- 
response shape as the filter clock frequency, number of terms in the 
Fourier series, and type of raised-cosine characteristic were varied. 

Croisier & Pierret [22] have shown theoretically that a digital 
modulator can be used to generate vestigial-sideband and phase- 
modulated signals which are identical to those obtained from analog 
techniques. They reported four types of digital modulators which 
have been produced, plus one experimental model. The production units 
generate either two- or four- level pulses with carrier frequencies of 
2.8 kHz to 64 kHz, while the experimental unit has eight-level pulses 
with a carrier frequency of 2.8 kHz. 


*The nonrecursive digital filter and the corresponding Fourier 
series expansion technique for its design will be described in Section 
2 . 1 . 
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Even though emphasis is being placed on the use of digital hard- 
ware in PAM systems, it can be seen from recent literature [23-26] 
that researchers are still investigating analog pulse-shaping network 
techniques . 


1.2.3 Related Developments 

Developments discussed previously were concerned primarily with 
the analysis and design of the PAM transmitter and receiver filters, 
since these are the areas of emphasis in this investigation. Two 
other areas which have become quite important in PAM system design 
are the use of encoding-decoding devices and adaptive equalization. 

Several authors [27-33] have discussed the application of coding 
techniques to increase the data transmission rate of binary data with- 
out a corresponding increase in the required bandwidth. Signaling at 
a rate higher than the theoretical Nyquist limit is accomplished by 
allowing controlled amounts of intersymbol interference. The encoding 
device for this method converts the original binary source sequence 
into a sequence of symbols {a^} which are applied to the transmitter 
filter. Encoding eliminates the interdependence of detected symbols 
at the receiver, so that a single detector error will not be propa- 
gated and result in an error burst. The method requires that the 
equivalent impulse response of the cascade network consisting of trans- 
mitter filter, channel, and receiver filter be of a special shape in 
order to control the allowable intersymbol interference. Kretzmer [31] 
has tabulated a number of these impulse responses which he calls 
partial- response pulses. By allowing intersymbol interference in the 
system impulse response, more than two possible sample levels exist at 
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the receiver. Consequently, a decoding device is used to generate 
binary estimates from the multilevel samples. Since more than two 
levels must be detected at the receiver, this technique is more 
sensitive to noise. • Nevertheless, the possibility of an increased 
data rate without increased bandwidth is a definite advantage for 
certain low-noise environments. 

Lucky, et al . [34-47] have developed the technique of adaptive 
equalization for use with time-varying channels such as are found in 
the switched-line telephone network. High-speed data transmission 
over channels with only preset equalizers at the transmitter and/or 
receiver cannot be sustained for widely time-varying channel char- 
acteristics. Consequently, an automatic equalizer, consisting of an 
adjustable transversal filter and associated control circuitry, has 
been developed to modify automatically the filter tap coefficients. 
Various algorithms exist for determining the proper tap settings. 

One technique is to transmit a test signal through the cascaded 
channel and automatic equalizer during a so-called "training period".. 
Typical test signals include a sum of sinusoids equally spaced 
throughout the channel bandwidth, a sequence of isolated pulses, and 
a wide-band pseudo-noise signal. The equalizer output response is 
compared to a locally generated reference response and the corre- 
sponding error signal used to increment the tap settings so as to 
reduce error. A more recent technique [37], called the "decision- 
directed method", is to continually adjust the tap settings during 
the transmission of actual data. In this scheme the detected symbols 
at the receiver are assumed to correctly represent the input symbols 
(which should be a valid assumption if the equalizer is in 
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near-optimum adjustment) thus eliminating the. need for a locally- 
generated reference. The cross correlation between the regenerated 
input symbols and the actual received waveform is used to determine 
the error signal for tap-setting purposes. Problems exist in the 
initialization of either of these tap-setting algorithms. The test- 
signal technique operates on the principle of peak-distortion mini- 
mization and the algorithm will not converge if the initial peak 
distortion is greater than some critical value. The decision- 
directed technique, which minimizes mean-square distortion, takes 
an abnormally long time to converge if the initial tap settings do 
not place the equalizer in near-optimum adjustment. 

1.3 OUTLINE OF STUDY 

Problems associated with PAM baseband data transmission were 
reviewed in Chapter 1. The remaining chapters are devoted to the 
design of pulse-shaping digital filters to solve some of these 
problems. 

Chapter 2 concentrates on the class of digital filters to be 
employed as pulse-shaping devices, i.e., the finite-duration impulse 
response filter. This class, which consists of the digital trans- 
versal filter and the frequency-sampling digital filter, can be 
implemented in hardware or by a general-purpose computer algorithm. 
Block diagrams, that would be used for a hard-wired digital- logic 
implementation, and difference equations, that would be programmed 
for a general-purpose computer implementation, are included for each 
filter. The design techniques reported in the literature, which are 
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mainly frequency-domain procedures, are reviewed and design examples 
are included. 

An optimization algorithm is developed in Chapter 3 to design 
the frequency-sampling digital filter from impulse-response specifi- 
cations. This technique is applicable to design digital filters to 
shape pulses for data transmission over unknown channels or .for 
matched-filter detection of pulses in noise. Equations are derived 
for the filter impulse response as a linear function of unknown 
multiplier coefficients, which must be determined to complete the 
filter design. It is then shown how these equations can be employed 
to constrain the filter impulse response to meet certain time-domain 
specifications. A linear-programming optimization algorithm is used 
to solve the constraint equations for the unknown filter coefficients. 
Details regarding linear programming and its computer implementation 
for digital filter design are found in the appendices. The procedures 
for designing digital PAM transmitter filters and digital matched 
filters are outlined by the use of examples. 

Application of the optimization algorithm to the design of 
digital filters for equalization of data channels that have known 
characteristics is developed in Chapter 4. Equations which are 
linear functions of the equalizer multiplier-coefficients are derived 
for the equivalent impulse response of the cascaded digital-filter 
equalizer, D/A converter, channel, and receiver filter combination. 
These equations are then used in the manner of Chapter 3 to constrain 
the received pulse shape so that the linear-programming algorithm can 
be used to specify the digital equalizer coefficients. Examples are 
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included for the design of both transversal filter and frequency- 
sampling filter equalizers. 

Finally, in Chapter 5 conclusions are drawn from the results 
of this investigation and problems for future study are recommended. 



CHAPTER 2 


FINITE-DURATION IMPULSE-RESPONSE FILTERS 

A class of digital filters characterized by a finite number of 
nonzero values in the impulse-response will be considered in this 
chapter. This class of filters has particular advantages for time- 
domain design techniques since the impulse response need be specified 
for only a finite time interval. Specifically, the impulse-response, 
which is denoted as h(mT), is defined as a set of N values 
{h„ ,h, , . . . ,h„ .. } with h = 0 for all other values of m. Two methods 
exist for the realization of this filter type, viz., the digital 
transversal filter and the frequency-sampling digital filter, which 
are nonrecursive and recursive implementations respectively. 

2.1 DIGITAL TRANSVERSAL FILTERS 

The digital transversal filter is a straightforward digital 
implementation of the analog transversal filter that was originally 
introduced by Kallman [48]. The use of digital logic affords more 
precise filtering at a lower cost and smaller size [19,49] than is 
possible with the tapped delay line and resistive summing network 
of the conventional analog transversal filter. 

2.1.1 Transversal Filter Characteristics 

A block diagram of the digital transversal filter is illustrated 
in Fig. 2-1. Each square signifies a delay element of T seconds, so 
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a collective delay of NT seconds can be realized by an N-stage shift 

register. The multiplication and summation elements are constructed 

from digital logic circuits. It should be noted that each of the 

signal samples, x(nT) and y(nT), is actually an M-bit digital word, 

M 

where 2 is the number of quantization levels used in the digital 
representation. Consequently, each set of time delays is realized 
by a parallel bank of M shift registers and the arithmetic elements 
multiply or add M-bit digital words. It can be observed from 
inspection of Fig. 2-1 that the difference equation for the trans- 
versal filter is written as 

N-l 

y (nT) = ) h x(nT - mT) . (2.1) 

L i m 

m=0 


Taking the z-transform of (2.1) gives the corresponding transfer 
function 


Y(z) A 
X(z) 


H(z) = 


I 

m=0 


h 

m 


-m 

z 


( 2 . 2 ) 


The sequence 


x(nT) 


1 

0 


n = 0 
n ^ 0 , 


(2.3) 


applied to a digital filter corresponds to the unit impulse which is 
applied to a continuous filter, because the response to this particu- 
lar sequence is the inverse z-transform of the filter transfer func- 
tion. Hence, the reason for setting the transversal-filter multiplier 
coefficients equal to the desired unit-impulse response hQ,...,h 
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can be seen by applying the sequence {1,0,0,...} to the network of 

Fig. 2-1. It will be assumed that there are no delay elements prior 

to the first multiplier tap. Hence, the first of the N impulse- 

response values will occur at time zero, i.e., h = h(mT). It will 

m 

also be assumed that x(n-T) = 0 for n < 0. 

Since (2.1) is the discrete equivalent of the convolution inte- 
gral as applied to continuous linear systems, the nonrecursive reali- 
zation is sometimes called a "direct-convolution filter". Stockham 
[50] has shown that by using the f ast-Fourier-transform (FFT) algo- 
rithm to calculate y(nT), the required computation time can be 
reduced significantly for sufficiently large N. When the FFT is 
used in this manner, the realization is called a "fast-convolution 
filter". 

In addition to having a finite impulse response, other dis- 
tinctive characteristics of the transversal filter include the 
absence of stability problems in the realization, a piecewise- 
linear phase characteristic which generates a uniform time delay 
and the need for a large number of coefficient multipliers to obtain 
sharp amplitude characteristic cutoffs. Since the transfer function 
(2.2) has only zeros there is no possibility of poles falling outside 
the z-plane unit circle; hence, the filter is inherently stable. The 
amplitude and phase characteristics will be illustrated shortly. 

2.1.2 Fourier-Series Design Technique 

Several conventional methods for the design of nonrecursive 
filters have been reviewed by Kaiser [51]. These techniques, which 
are all frequency-domain procedures, have in the past formed the 
basis for transversal filter design. Perhaps the most useful has 
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been the Fourier series expansion technique. Since the desired 
digital filter frequency characteristic, H(f), must always be 
periodic, it can be expanded in a Fourier series over the range 
| f | < 1/2T , i.e. , 


H(f) = -zr~ + V [a cos(2irnTf) + b sin(2imTf ) ] 
/ L i n n 

n=l 


where 


and 



(2.4a) 


(2.4b) 


(2.4c) 


The Fourier-series coefficients can be shown to represent the desired 
impulse-response values. Consequently, the desired filter coeffi- 
cients are obtained by: (1) replacing the sinusoidal functions with 

their Euler expansion equivalents, (2) making the standard substitu- 
tion 


z 


j 2frfT 


(2.5) 


and (3) comparing the resulting transfer function with (2.2). Of 
course, the series must be truncated to obtain a finite number of 
coefficients, and in general the filter will have complex-valued 
coefficients and will require future values of the input sequence, 
i.e., the summation in (2.2) will include negative values of m. 

As an example, Kaiser has considered the problem of designing 
an ideal digital transversal differentiator. The Fourier-series 
coefficients which correspond to the desired frequency characteristic 
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H(f) = J 2 Trf |f j £ ~ , (2.6) 

can be found by substituting (2.6) into (2.4b) and (2.4c). giving 

a=0 n = 0 , 1, . . . ,N 

n 

and (2.7) 

b n = (-l) n+1 n = 1,2,..., N . 

n nl 

The transfer function is obtained by Substituting (2.7) into the 

exponential form of (2.4a) with the series truncated to N =' 9. 

-9 

Then the transfer function is multiplied by z to eliminate the 
need for future input values which gives a transfer function with 
18 non-zero coefficients, i.e., 

H(z) = E ~ nr [zn ~ 9 ~ z_n ~ 9] • (2 - 8) 

n=l 

Characteristics of the transversal filter corresponding to (2.8) 

with T = Is are given in Fig. 2-2.* The Gibbs phenomenon evident 

in the amplitude characteristic is due to the series truncation, 

which illustrates the transversal-filter property that a large 

number of taps are required for sharp amplitude cutoffs, i.e., more 

terms are needed in the series expansion. As expected, the phase 

characteristic is exactly linear over the entire bandwidth. The 

linear phase shift which has been subtracted from the ideal 90° 

-9 

phase shift results from the z term in (2.8). It can be seen 


*A11 illustrations of digital filter frequency- and time-domain 
characteristics were generated with the ZFRAPP and ZTRAMP subroutines 
which are described in the technical report by Routs & Burlage [52]. 
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(a) Amplitude Characteristic 



(b) Phase Characteristic 



(c) Impulse Response 
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Frequency- and Time-Domain Characteristics of Digital 
Differentiator Designed by Fourier Series Technique. 
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from inspection of Fig. 2-2 (c) that the impulse response contains 
the same 18 values predicted by (2.8) plus the zero value at 9T 
corresponding to a^ = 0. 

2.1.3 Window- Function Design Technique 

Reduction of the Gibbs phenomenom has been studied by Kaiser 
and others [53-56] resulting in the development of the so-called 
"window function" design technique for transversal filters. First, 
a Fourier series expansion of the desired frequency characteristic 
is obtained. This series is then multiplied by the window function, 
which is time limited and is used to modify the series coefficients 
so as to reduce the normal truncation error. Since multiplication 
in the time domain is equivalent to convolution in the frequency 
domain, multiplication of the time-domain coefficients by a proper 
window function results in a smoothing of the sharp transitions 
found in the desired frequency characteristic, which then can be 
more closely approximated. The modified Bessel window function 
discovered by Kaiser, i.e., 



is shown in Fig. 2-3 (a) and has been used to redesign the ideal 
differentiator. Multiplying each coefficient of (2.7) by a cor- 
responding window-function value w(nT), with t = 9T, results in 
a transversal filter which has the characteristics shown in Fig. 
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2-3. It is obvious that the ripple effect, so evident in Fig. 

2-2 (a), has been reduced significantly. Discussion of the advan- 
tages, disadvantages, and tradeoff possibilities of various window 
functions are found in the previously cited references.. 

Frequency-domain design techniques developed subsequent to 
Kaiser [51] have been primarily applications of optimization algo- 
rithms to minimize errors between the desired and actual filter 
frequency characteristic. Since these techniques are applicable 
to both the transversal and frequency-sampling filter realizations, 
discussion will be deferred until the frequency-sampling filter has 
been described. 

2.1.4 Time-Domain Design Techniques 

Techniques for designing digital transversal filters in the 
time domain have not received much emphasis in recent years nor have 
they been developed to the degree of the frequency-domain methods. 
The majority of the work that has been done is mainly applications 
of classical numerical analysis procedures for prediction, smoothing, 
and differentiation. Monroe [57], in summarizing these procedures, 
has derived equations to define the N coefficients, h , of (2.1) so 
that the designed filter has a maximum noise reduction capability 
while simultaneously having desired differentiating and/or pre- 
dicting characteristics. The input to the filter is assumed to be 
a polynomial of finite degree with additive white noise. He has 
also developed a method to specify the coefficients {h^} so that 
the filter step response will meet desired criteria, such as over- 
shoot, and still have noise rejection capabilities. 
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Cavin, Ray, and Rhyne [58] have considered the time-domain 
design of an optimal transversal filter which would give an impulse 
output when an expected seismic waveform, called a Ricker wavelet, 
is applied to its input. They studied three error criteria for 
indicating the difference between the actual filter response and 
the ideal impulse, namely, a weighted-mean-square-error, a weighted- 
absolute-error, and a minimax error. Filter coefficients that 
would minimize each of these error functions were determined with 
a linear programming algorithm. One problem with this technique 
was that all of the resulting designs were of the high-pass fre- 
quency characteristic, and could thus amplify noise components 
that might be present in the seismic waveform. 

2.2 FREQUENCY-SAMPLING DIGITAL FILTERS 

The frequency-sampling filter, a recursive realization of 
the finite-duration impulse response filter, was first introduced 
by Rader & Gold [59]. The name of the filter results from their 
frequency-domain design technique rather than an implied descrip-, 
tion of the filter realization or operation. 

2.2.1 Frequency-Sampling Filter Characteristics 

A block diagram of a frequency-sampling filter is illustrated 
in Fig. 2-4. The network consists of a comb filter in series with 
a parallel combination of 1 + [N/2]* digital resonators each 


* The notation I [N/2] will be used to indicate the integer 
portion of N/2. 
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Fig. 2-4. Frequency-Sampling Digital Filter. 
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"tuned" to a different frequency. The reasons for this terminology 
will become evident when the filter frequency characteristics are 
considered. First, however, the filter difference equations and 
transfer functions must be derived. 

The frequency-sampling digital filter shown in Fig. 2-4 is 
described by the following difference equations: 

m(nT) = x(nT) - x(nT - NT) , 

H 0 

y Q (nT) = — m(nT) + y Q (nT - T) , 

1c 

y k (nT) = - 21^ [m(nT) - cos(2irk/N) m(nT - T) ] 

- 2 cos (2irk/N) Y k (nT - T) + Y k (nT - 2T) 

k = 1,2,. . . , I [N/2 ] , (2.10c) 

and 

I [N/2] 

y (nT) = y k (nT) • (2 . lOd) 

k=0 

An overall transfer function could be derived by taking the 
z-transform of each of these equations, as was done for the trans- 
versal filter in (2.2), and then combining the resulting expres- 
sions. Alternatively, the approach of Gold & Jordan [60] will 
be followed in order to better illustrate certain properties of 
the filter. The discrete Fourier transform (DFT) of the finite 
impulse response sequence {h^ , . . . ,h^_^} yields 


(2.10a) 

(2.10b) 
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N-l 

H(k«) = y h e _;3mTkQ k = 0,1,..., N-l , (2.11) 

/ . m 

m=0 

where = 2 tt/NT is the increment between samples in the frequency 
domain. Those samples, H(kft), correspond to evaluation of H(z). at 
N equally-spaced points around the unit circle in the z-plane, and 
consequently are values of the continuous frequency response for 
the filter. Taking the inverse DFT of these frequency samples 
yields an expression for the impulse response as a function of the 
samples, namely 

N-l 

h m = ^ ^ H(kO) e jkfimT . (2.12) 

k=0 

Hence, by substituting (2.12) into (2.2), interchanging summations, 
and using the closed-form expression for the geometric series, the 
following transfer function is obtained which is also a function of 
the frequency samples. 


H(z) 



I 


H(kfl) 


l-eJ ™.' 1 


(2.13) 


This function could be realized with a network very similar to 

-N 

that of Fig. 2-4. The term 1 - z represents the transfer func- 
tion for the comb filter and the terms H(kfi)/N, which may be 
complex valued, represent the output multiplier coefficients for 
the resonators. The primary difference is that each resonator 
would be of first-order with a single complex-coefficient multi- 


plier. By assuming that 
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H(kfl) = H*(Nf2 - kfi) , . (2.14) 

complex conjugate terms in (2.13) can be combined in the manner 


1 - 


H(kfl) 

jTkfl -1 
e z 


H(NQ ~ kO) 

, jT(N-k)n -1 
1 - e z 


(2.15) 


2 Re{H(kft)} - 2[Re{H(kfl)> cos(Tkfl) + Im{H(kfl) }■ sin(Tkfl) ] z 

1-2 cos (TkQ) z ^ + z ^ 


For real H(kfi) 
to Fig. 2-4, i, 


H(z) = 


1 - . z 
N 


-N 


(2.15) becomes the transfer function corresponding 
e. , 



+ 


I [N/2 ] 

I 


k=l 


(-1)^211 [1 - cos(2mk/N)z 
1-2 cos (2irk/N) z ^ + z ^ 


(2.16) 


where real H(kn) are denoted by H^. The real frequency samples 
have been multiplied by (-1) since Rader & Gold [59] have shown 
that this multiplication is necessary to have a smooth interpola- 
tion through the frequency samples by the continuous frequency 
response of the filter. Although the transfer function (2.16) 
holds for all odd N, it is only valid for even N if = 0. If 

this is not the case, the function can be corrected by replacing 
the term the summation with (-1)^^ / (1 + z ^) . 

2.2.2 Conventional Frequency-Domain Design Techniques 

All techniques that have been used to design frequency- 
sampling filters have been frequency- domain methods. The original 

technique of Rader & Gold [59] was to specify the filter multiplier 
Ic 

coefficients (-1) H^N directly from the samples of the desired 
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continuous amplitude frequency characteristic taken at N equally- 
spaced frequencies. An insight into this design technique can be 
obtained by inspection of Figs. 2-5 and 2-6. Significant poles 
and zeros of (2.16) are illustrated in the z-plane of Fig. 2-5 
for a three-resonator filter with N =20. The comb filter results 
in 20 equally spaced zeros around the unit circle.. Therefore, as 
frequency is increased from zero to the half-sampling frequency, 

1/2T, corresponding to movement midway around the unit circle,. the 
comb-filter magnitude will appear as a series of 10 equal-magnitude 
lobes or "teeth" with the zero value between adjacent lobes resulting 
from a z-plane zero. On this semicircle, each resonator is repre- 
sented by a complex pole which cancels the effect of the comb-filter 
zero at that resonant frequency. Consequently, those frequency 
components contained within a bandwidth centered, on the resonant 
frequency are passed by that resonator stage, and the overall filter 
frequency characteristic is determined by the sum of these individ- 
ual passbands. 

A typical set of frequency samples which would correspond to 
the pole-zero plot of Fig. 2-5 are given in Fig. 2-6. This design 
of a filter with three resonator stages and a comb filter with 20 
stages is an attempt to obtain an ideal low-pass frequency charac- 
teristic. Frequency- and time-domain characteristics for the fil- 
ter are plotted in Fig. 2-7. It is evident from the amplitude 
characteristic that the continuous response passes through each of 
the specified frequency samples, even though the ideal character- 
istic could not be obtained. The piecewise linear phase shift 
over the entire frequency range and the finite impulse response 
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are distinctive characteristics of this type of filter as they were 
for the transversal filter. 

The reason for the impulse response being finite can be seen by 
applying the input sequence {1,0,0...} to the network of Fig. 2.4. 

The comb filter causes two sequences to be applied to each resonator, 
i.e., the original sequence and a delayed negative version NT seconds 
later. Then by taking the inverse z-transform of a resonator trans- 
fer function to obtain 

. 2H 

H (mT) = (-1) cos(2TTkm/N) , (2.17) 

K IN 

it is obvious that the resonator impulse response repeats at NT 
seconds. Hence, the delayed negative impulse cancels the effect of 
the original impulse giving a zero output for all times after NT 
seconds. 

Any implemented digital filter requires multiplier coefficients 
to be represented by digital words having a finite number of bits. 

This coefficient quantization error can have a serious effect for the 
frequency-sampling filter which requires poles to be placed directly 
on the unit-circle stability boundary. Weinstein [61] and others who 
have studied this problem have shown that acceptable performances can 
be obtained by moving all the critical poles and zeros slightly inside 
the unit circle, i.e., at a radius of r = 1 - 6. The value for 6 can 
usually be in the range of 2 ^ to 2 with little noticeable change 
in the filter characteristics. 

2.2.3 Optimization Technique for Frequency-Domain Design 

Recently developed procedures to design frequency-sampling fil- 
ters [62-64] employ optimization algorithms to minimize the maximum 
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frequency response deviation of the designed filter from the ideal 
response over a specified frequency range. The basic function used 


in the optimization procedures can be derived by making the substitu- 
tion z = e^ W ^ in (2.13) and manipulating the resulting expression to 


yield 


H(u>) 


1 -j [ (uNT/2) (1 
N 


N-l 

1/N)] y 
k=0 


H k e-jf Tk /Nsin( a)NT /2) 
sin(a)T/2 - mk/N) 


(2.18) 


This is an equation for the continuous frequency response of the fil- 
ter as a linear function of its real frequency samples H^. Hence, it 
is possible to apply linear optimization techniques to select a set 
of values to minimize or maximize some characteristic of the con- 
tinuous frequency response. Gold & Jordan [62] have used this tech- 
nique to minimize the maximum sidelobe amplitude in the out-of-band 
response for a low-pass filter. A frequency characteristic resulting 
from their design procedure is compared in Fig. 2-8 with a character- 
istic resulting from the original procedure. Both filters have 
N = 256, T = 1 s, 1^ = 1 for 0 < k < 31 and = 0 for 35 < k < 128. 
The frequency samples ^ 33 ’ anc ^ ^34 an (2.18), which were set to 

zero for the original design, define a transition band whose values 
are optimally selected in the second design as 0.7, 0.225, and 0.02, 
respectively. By using the optimization technique, the maximum out- 
of-band response is reduced from -28 dB to -85 dB. 

Both the original and optimum frequency-sampling design proce- 
dures can be used to design transversal filters. After obtaining 
the set of H^'s by either method, the inverse DFT, (2.12), is used 
to calculate the corresponding transversal filter multiplier 
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Fig. 2-8. Improvement in Stop-Band Aaplitude Response of Low-Pass 
Flltar Due to Optimization Technique. 
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coefficients {h^}. Unfortunately, this approach results in a trans- 
versal filter with N multipliers, 256 for the above example, and 
consequently a more inefficient design is obtained than with the 
recursive realization. This illustrates one of the main advantages 
of the frequency-sampling filter over the transversal filter, i.e., 
fewer multipliers are needed to obtain sharp cutoff-frequency charac 
teristics. 

No technique to go directly from time-domain specifications to 
filter coefficients has been reported for the time-domain design of 
the frequency-sampling filter. This type of time-domain design is 
the subject of the next chapter. 



CHAPTER 3 


DIGITAL FILTER DESIGN FROM IMPULSE-RESPONSE SPECIFICATIONS 

Many of the results which have been derived for optimum PAM 
data transmission, as outlined in Section 1.2.1, were expressed as 
equations for the required impulse response of the transmitter and 
receiver filters. However, no straightforward technique has been 
developed to design and synthesize filters directly with digital 
hardware from these time-domain specifications. An algorithm is 
introduced in this chapter for that purpose. The procedure consists 
of three steps. First, equations are derived for the impulse response 
of a digital filter as a linear function of its unknown multiplier 
coefficients. Next, the impulse response is constrained with these 
equations to have appropriate values at certain critical time points. 
Finally, a linear programming routine is applied to solve the con- 
straint equations and thus define the unknown filter coefficients. 

3.1 TIME-DOMAIN DESIGN PROCEDURE 

Expressions are derived in this section for the frequency-sampling 

filter impulse response, h , as a function of its unknown coefficients, 

m 

H^, k = 1,2 , . . . ,1 [N/2] . This derivation is not needed for the trans- 
versal filter since its multiplier coefficients are equal to the speci- 
fied impulse-response values. The inverse DFT of a sampled frequency 
characteristic, which was obtained in Section 2.2, is an equation for 
the frequency-sampling filter impulse response, namely 
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N-l 

h m = | ^ \ e jkfimT 111=0,1 N-l. (3.1) 

k=0 

A z transfer function (2.16) was derived by assuming that the fre- 
quency samples were real and possessed the symmetry 
These restrictions were assumed to simplify the design problem. If 
the were not real an optimization algorithm with complex arith- 
metic would be required. Furthermore, if the symmetry did 

not hold, filters with complex impulse responses would be possible. 
With these two assumptions (3.1) can be simplified. However, the 
cases for odd and even values of N will be considered separately. 


\ “ Vie’ 


3.1.1 Odd Number of Impulse-Response Values 

Typical frequency samples for the N-odd case are illustrated 
in Fig. 3-l(a). If the H Q sample is ignored, the samples are sym- 
metrical about the line k = N/2. Therefore (3.1) can be written as 

„ C-D/2 „ 

^=-0+ £ _k [e j2»km/N + e j2,(N-k)m/Nj _ (3 . 2 

k=l 

and then simplified to give 


h 

m 


N 


+ 


(N-l)/2 

I 


k=l 


2H k 

— — cos(2iTkm/N) 


(3.3) 


which is a real, linear function of the multiplier coefficients. 

It can be seen from (3.3) that the impulse response for this 

case has the symmetry property h^ = h^_ m for m = 1,2 , . . . , (N-l) /2 , 

1c 

If the multiplier coefficients are multiplied by (-1) , as was done 
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in (2.16) to obtain a smooth interpolated frequency response, the 
resulting impulse response 


h 

m 


H, 


N 


+ 


(N-l)/2 

i ( - i>k 

k=l 


2H k ' 

— — cos (2iTkm/N) . 


(3.4) 


still has the symmetry h^ = h^ . However, the impulse- response peak 
for (3.3) occurs at m = 0, whereas for (3.4) the peak would occur at 
m = (N±l)/2. This is evident by comparing Figs. 3-2 (a) and 3-2 (b), 
where the individual resonator impulse-response envelopes for two 
filters are shown. For the first filter, which corresponds to Fig. 
3-l(a) and does not contain the (-1) factor, the individual enve-r 
lopes are all in phase at t = 0. Conversely, the envelopes for the 
second filter which has the (-1) coefficients are in phase at 

t = NT/2. 

3.1.2 Even Number of Impulse-Response Values 

Frequency samples for the N-even case, as illustrated in Fig. 
3-l(b), have an axis of symmetry at N/2 which is also the location 
of a frequency sample. Consequently, both the k=0 and k=N/2 samples 
must be considered separately in writing the impulse response equa- 
tion. With these observations, (3.1) becomes 


h 

m 


H r 


N 


+ 


(N/2)-l 

V k 2H k 

y (-1) — rc cos (2ukm/N) + 

k=l 



cos (nnr) 


V2 

N 


(3.5) 


The inclusion of the (-1) term again results in a shift of the 
impulse-response peak; this time from m = 0 to m = N/2. The impulse 
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response symmetry = h^_ m also holds, as can be seen by inspecting 
Fig. 2-7(c) where the peak occurs at N/2 = 10. 

The signal component illustration of Fig. 3-2 also provides 
insight into the problem of designing frequency-sampling filters to 
have specified impulse responses. In a time-domain analysis, multi- 
plier coefficients {H^} correspond to the components present in a 
Fourier series expansion of the impulse response. They are the weight 
of the sinusoidal Fourier components plotted in Fig. 3-2. Hence, the 
linear-programming routine can be thought of as a way to determine 
which Fourier components are required and what their corresponding 
weights should be so that the sum of the resulting sinusoids will 
give the desired impulse response. Of course, these component values 
are evaluated and added only every T seconds resulting in the discrete 
filter output. 


3.1.3 Formulation of Constraint Equations 

Equations (3.4) and (3.5) are discrete-valued expressions which 
can be used for the linear-programming time-domain design of the 
frequency-sampling filter. An expression for the unit-impulse 
response envelope, which is given by the sum of the individual reso- 
nator impulse- response envelopes, can be obtained from (3.4) or (3.5) 
by making the substitution t = mT in (3.4) and (3.5), i.e. 


h(t) 



(N-l)/2 

z 


k=l 


(-D 


k^k 

. N 


cos (2irkt/NT) 


(3.6) 


and 



-2 


V 



Fig, 3-2. Comparison of Impulse-Response Components for Two 
Sets of Multiplier Coefficients. 
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H, 


(N/2)-l 2R 

i(t) = ^ + ("D k -5^ cos (2-rrkt/NT) + (-1) N/2 cos(irt/T) 


N 


N 


k=l 


(3.7) 


These expressions are continuous over the required range 0 <_ t _< NT 
and are valid for N odd and even respectively. The envelope is a 
good approximation to the continuous signal that would be generated 
by a D/A converter and smoothing filter at the digital filter output. 
Since continuous linear functions of the multiplier coefficients are 
now available, the impulse-response slope and other higher deriva- 
tives of h(t) can be constrained. For example, the required slope 
equation for N odd is 

(N-l)/2 

h' (t) = y k(-l)k + "*' h sin (2iTkt/NT) . (3.8) 

NT ^- J k 

k=l 


Equations (3.4) through (3.8) have been written in a form such 
that the theory of linear programming can be applied to solve for 
the unknown H^, k = 0,l,...,I[N/2]. A typical set of equations which 
might be used to specify a desired impulse response can be stated as 
follows. The slope of the impulse response at time t^, i.e., 


h'(t 1 ) = 


11 0 


H„ + 


*12 H 1 


+ a lR H R-1 


(3.9a) 


is to be minimized subject to the constraint equations 


h(t l> - a 21 H 0 + a 22 H 1 + ••• + a 2R H R-1 i 1 

and 

h(t 2 ) = a 31 H q + a 32 4- ... + a^ R H R _^ = 0 • (3.9c) 
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The coefficients a„ are knoxm, e.g., a 22 = (-2/N) cos (2irt^/NT) . The 
for k = R to I[N/2] are set to zero. Emphasis in this study will 
be on the design of low-pass and band-pass filters for which, 
frequency samples in certain regions are zero. Hence, the parameter 
R is defined where the frequency samples from k = R to k = I[N/2] 
are set to zero during the constraint equation formulation. The use 
of linear programming to solve for the set of R unknown multiplier 
values {H^} will be discussed in the two remaining sections of this 
chapter. 

3.2 LINEAR-PROGRAMMING OPTIMIZATION ALGORITHM 


Linear programming is an iterative procedure to solve for the 
set of R non-negative variables {x^} which maximizes or minimizes 
the linear function 


Z = c^ + c 2 x 2 + . . . + c R x R , 


(3.10) 


subject to the M linear constraint equations 


a il X l + a i2 X 2 + ' • ' + a iR x R 


i = I b . i = 1,2 


M 


(3.11) 


Without loss of generality, it will be assumed that all b^ >_ 0 and 
that Z is only to be maximized. Equation (3.10) is called the 
"objective function". A solution to (3.11) for which all x^ . >_■ 0 is 
called a "feasible solution". A feasible solution with no more than 
M positive x. is called a "basic feasible solution", and a basic 

l 

feasible solution which maximizes (3.10) is called an "optimal basic 
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feasible solution". A solution for which Z can be made arbitrarily 
large is called an "unbounded solution". 

The systematic procedure to solve (3.10) and (3.11) for the 
optimal basic feasible solution, called the "simplex method", was 
first developed by Dantzig [65] and has been discussed in detail 
by Gass [66] and Hadley [67]. The simplex method is a procedure 
which consists of starting with a basic feasible solution to (3. 11) 
and then moving from the first basic feasible solution to a second, 
etc. each time moving to a new solution which increases the value 
of Z. After a finite number of steps, an optimal basic feasible 
solution will be reached, i. e. , no other solution yields an increase 
in Z, The basic mathematical technique of the simplex procedure and 
a computer implementation of the procedure, called LINPO, are pre- 
sented in Appendices A and B, respectively. Two solutions to the 
following example problem are also included in the appendices to 
illustrate the simplex method. A hand-computation solution is 
developed in Section A. 2, while a LINPO subroutine solution is pre- 
sented in Section B.2. 

Consider the problem of finding a basic feasible solution which 
will maximize the objective function 

Z = + x 2 (3.12) 

subject to the constraints 


X, 

+ 

x 0 > 2 

1 


2 — 

3x 

+ 

x„ < 15 

1 


2 — 

x. 

+ 

2x 2 1 9 

1 




and 
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Two or three variable linear-programming problems can be repre- 
sented and solved graphically. Since a feasible solution requires 
positive variables, all possible solutions to the problem must be 
in the first quadrant of the x^-x^ coordinate system as shown in 
Fig. 3-3. The three constraint equations further restrict the 
solution to the interior of the crosshatched region. Consequently, 
all points inside the region represent feasible solutions. In order 
to determine the optimal solution, a family of straight lines has 
been superimposed on the figure with each line representing points 
that satisfy (3.12) for one value of Z. Hence, the optimal solution 
is that point or set of points located in the feasible region which 
also lie on the objective function line with the largest value of Z. 
For this example, the optimal solution is denoted by the point (4.2, 

2.4) where Z = 6.6. 
max 

3.3 DESIGN EXAMPLES 

Two examples are presented to illustrate the applicability of 
the LINPO algorithm to the time-domain design of frequency-sampling 
filters. One example involves the design of a pulse-shaping digital 
filter to generate a raised-cosine pulse and the second involves the 
design of a partial-response pulse-shaping filter. Although the two 
design examples involve pulse-shaping filters, both can be construed 
to be designs of digital matched filters for the detection of the 
corresponding pulse in a background of additive white Gaussian noise. 

The general procedure to be used for LINPO designs from impulse 
response specifications is outlined in the following table. 
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TABLE 3.1 

LINPO ALGORITHM DESIGN STEPS 


1) Determine the critical time points of the desired impulse 
response, write constraint equations for each point, and deter- 
mine an objective function. Since finite impulse-response filters 
are being designed, only a finite number of points of an infinite 
response specification can be constrained. Consequently, the 
most significant portion of the response must be ascertained and 
only points in that time interval constrained. For PAM waveforms , 
this interval is centered on the pulse peak and extends in both 
directions for a few Nyquist intervals. 

2) Determine N, the number of impulse- response values, and T, the 
filter sampling interval. Time-domain restrictions on these 
values are that NT must equal the time interval corresponding 
to the significant portion of the impulse response and NT/2 
must equal the time of the pulse peak. Physical restrictions 
are that the digital-logic speed limits the minimum value of T 
and the feasible number of shift-register stages for the comb 
filter limits the maximum value of N. 

3) Determine R, the maximum number of resonator stages in the filter. 
If the bandwidth of the desired impulse response is known, R is 
defined since R/NT equals that bandwidth for frequency-sampling 
filters. Otherwise, the bandwidth is estimated from the impulse 
response shape and that estimate is used to calculate an initial 
R. Restrictions are that R <_ I[N/2] and that the amount of dig- 
ital logic required for the resonator stages is feasible. 

4) Calculate coefficients of the unknown multiplier values in the 
constraint equations and objective function for input to the 
LINPO subroutine. Solve for the set { } with the LINPO sub- 
routine. 

5) Calculate the impulse response of the filter defined by the 
{H^} from the LINPO design. If the response is unsatisfactory, 
change the constraint equations and/or the parameters N, T, and 
R so that the response is improved when the previous step is 
repeated. 


3.3.1 Raised-Cosine Pulse Filter 

A pulse shape frequently used for baseband data transmission 
is the raised-cosine pulse [6] which is described by 
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g(t) 


sin(nt/t. ) cos (cnrt/t, ) 
b b 

’ C/t b 1 - (2«/t. ) 2 ’ 

b 


(3.14) 


where t, is the bit time and 
b 


g(t) has the non-zero frequency spectrum 


G(f) = 


4- 0 < |f| < (l-a)f 

1 


(3.15) 

(1-a) f 1 < | f | 1 (l+a)f 1 . 


The nominal cutoff frequency f^ is equal to l/(2t^). A delayed ver- 
sion of (3.14), g(t - 3t^) is plotted in Fig. 3-4 for a =1. Only 
the portion of (3.14) in the interval of ±3^ on either side of the 
pulse peak has been considered significant since beyond this time 
interval the maximum amplitude is less than 0.002. The primary 
advantage of this pulse shape, as was discussed in Section 1.1, is 
the presence of the zero crossings at multiples of t^. Zero ampli- 
tude at these points eliminates the intersymbol- interference trans- 
mission impairment. The design problem therefore is to specify a 
filter which has its impulse response given by (3.14) 

The frequency-sampling filter impulse-response envelope, (3.6) 
or (3.7), is symmetrical about the point NT/2; so the response need 
only be constrained on one side of the peak. Consequently, a set 
of constraint equations which might be sufficient for the design are 


h (3t, ) = 1 h(5t, ) = 0 h'(4t,) < 0 

b b b — 

h(4t, ) = 0 h(6t, ) = 0 h'(5t.) < 0 

b b b — 


(3.16) 




Filter Design (R-6). 
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where it is desired to maximize 

Z = h' (4t fa ) + h’ (5t fc ) . (3.17) 

This objective function along with the two slope constraint equations 
will attempt to force the slope to be zero at the indicated zero- 
crossing points. This modification of the raised-cosine pulse would 
be desired to reduce the effects of timing jitter which might exist 
in transmission systems. 

Parameters, other than the multiplier coefficients, which must 
also be specified are the number of impulse-response values, the 
sampling interval, and the number of resonator stages. The relation- 
ship of these parameters to the known and unknown multiplier coeffi- 
cients in the frequency domain is illustrated in Fig. 3-5 for N even. 
Since the spectrum of the raised-cosine pulse has a low-pass charac- 
teristic, the frequency samples from k = R to k = N/2 must be set to 
zero. Furthermore, the spacing between the samples, 1/NT, has been 
fixed by the requirement that NT/2 = 3t b , which is the pulse peak. 

It will be assumed for this example that the logic speed is suffi- 
cient for T = t, /10. Hence N = (3t,)(2/T) = 60, which is a realistic 
b d 

value for the number of shift-register stages in the comb filter. 

With a = 1 in (3.15), it can be seen that the raised-cosine spectrum 
extends to 2f^ = 1/t^. Consequently, it will be assumed that R/NT = 
2f^ which then yields R = NT/t^ =6. 

The filter design has now been reduced to a linear-programming 
problem with six constraint equations specified by (3.16), an objec- 
tive function given by (3.17), and six unknowns, {H^} . However, 
when these constraint equations are expressed in the form of (3.9), 
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the a„ coef ficients , i = 1,2,..., 7 and j = 1,2,..., 6 must be calcu- 
lated before linear programming can be used. An auxiliary subroutine, 
called C0CAL1, has been developed for this purpose and is described 
in Appendix C. In addition to calculating the a„ coefficients, this 
subroutine calculates and establishes all other arrays and parameters 
required for a linear-programming solution and then calls the LINPO 
subroutine which attempts to find an optimal basic feasible solution. 

The use of C0CAL1 and LINPO to solve (3.16) and (3.17) for the 
unknown is presented as an example in Section C.l. The results 
of that linear-programming solution and the corresponding impulse 
response of the digital filter are shown in Fig. 3-6(a). Although 
the pulse has the desired magnitudes at (3 ± n )t^, n = 0,1,2, 3, and 
the slope at the zero crossings has been forced to zero, it does not 
have the a = 1 raised^-cosine shape. Since zero crossings are missing 
at 4.5t^ anc j 5.5t^, two more constraints, i.e., 

h ( 4 . 5 1, ) =0 

D 

and 

h(5.5t. ) = 0 

D 

were added to the set of (3.16) and the problem solved again. The 
second LINPO design, as shown in Fig. 3-6(b), yields a digital filter 
with an accurate raised-cosine impulse response. The difference 
between the impulse response and the ideal pulse is less than 0.0012 
for each of the 60 response values, which is less than 0.12% of the 
peak pulse value. The sum of the slopes at the two zero-crossing 
points, 4t^ and 5t^, could not be driven to zero for the second 
design, as indicated by the optimum value of the objective function 




Fig. 3-6. Comparison of Raised-Cosine Filter Designs for Two 
Sets of Constraint Equations. 
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Z = -0.036. Consequently, for applications where jitter is a defi- 
nite problem the first design would be desired. Conversely, if 
jitter is negligible the second design will generate an excellent 
raised-cosine pulse and could be used. It may seem obvious that the 

zero crossings at 4.5t, and 5.5t, should have been constrained in 

b b 

the first place. However, it is desirable to keep the number of con- 
straint equations to a minimum which in effect reduces the hardware 
required. 

A situation which did not arise during this example is the case 
of no feasible solution existing for a specified set of consistent 
constraint equations. The number of equations was increased from 
six to eight during this design, but the number of unknowns was held 
to six. If the number of equations is continually increased in this 
manner without a corresponding increase in the number of unknowns , a 
point will be reached where the problem is overconstrained in that 
it does not have a feasible solution. Consequently, the number of 
unknowns would have to be increased to obtain a feasible solution. 
However, this increase requires more hardware in the filter since 
each unknown corresponds to a resonator stage. Hence, it may become 
necessary to make tradeoffs between hardware complexity and the 
number of constraint equations employed in the design. A set of 
constraint equations with no feasible solution is illustrated in 
the next example. 

3.3.2 Partial-Response Pulse Filter 

The raised cosine pulse shapes have zero magnitude at multiples 
of the time parameter t^ ; hence, zero intersymbol interference. 
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Another class of pulse shapes which permit intersymbol interference, 
but in a controlled amount, are the partial-response pulses [31] 
defined by , 


g(t) 


L 


k 

n 


sin [ir(t - nt b )/t b ] 
ir(t - nt b )/t b 


(3.18) 


The significant portion of a typical pulse is plotted in Fig. 3-7 (a) 
for the set of integer weights k_^ = = -1, k_^ = k^ = 0, and 

kg = 2. These weights result in a pulse which has intersymbol inter- 
ference at times ±2t b from the peak but not at other multiples of t b> 
The corresponding pulse spectrum is 


G(f ) = y ~ sin 2 (7rf/f 1 ) 0 1 l f l 1 • (3.19) 

It is evident from a comparison of (3.15) and (3.19) that by allowing 
this specified amount of interference, the total channel bandwidth 
required for the pulse is halved and the need for the channel to 
possess dc-frequency response is eliminated. 

A frequency-sampling filter which generates the partial- response 
pulse is designed in the same manner as the raised-cosine pulse fil- 
ter. With the peak of the impulse response specified by NT/2 = 5t, 

b 

and assuming that T = t b /5, the resulting number of impulse response 

values for the pulse is N = (5t,)(2/T) = 50. The number of resonators 

b 

stages, R, required to generate the pulse is obtained by setting 

R/NT = f^, which is specified in (3.19) to be the required bandwidth. 

It then follows that R = NT/2t, = 5. 

b 

Critical values of the partial-response pulse are the peak value 
of 2.0 at 5t b , the intersymbol interference magnitude of -1.0 at 7t b> 



h(mT) 



Fig. 3-7. Comparison of Partial-Response Filter Design to Ideal Pulse. 
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and the zero magnitudes at 6t^, 8t^ and 9t^. Although not a critical 
value, another distinctive characteristic of the pulse is the zero 
slope at 6.756t^. A corresponding set of constraint equations are 

h(5t, ) = 2 h(7t, ) = -1 h (9t, ) = 0 

b b b 

(3.20) 

h ( 6 1 , ) = 0 h (8t, ) = 0 h ' (6. 756t, ) < 0 

b b b — 

By maximizing the objective function 

Z = h'(6.756t, ) , (3.21) 

b 

subject to the last constraint equation, it should be possible to 
force the slope to zero at the desired time point. However, the 
linear-programming problem defined by (3.20) and (3.21) has six 
constraint equations but only five unknowns. As a consequence, 
application of the COCAL1 and LINPO subroutines to solve the prob- 
lem resulted in the statement "infeasible", i.e., no feasible solu- 
tion exists for five variables. Consequently, it was necessary to 
either increase the number of unknowns or decrease the number of 
equations. All of the constraints were desirable, so the number of 
unknowns had to be increased to seven before the satisfactory solu- 
tion shown in Fig. 3-7(b) was obtained. 

The difference between the designed filter impulse response and 
the partial-response pulse is less than 0.84% of the peak for each of 
the 50 response values. The price paid to obtain a feasible solu- 
tion, however, was an increase in pulse bandwidth from 5/NT to 7/NT 
and the possible addition of two resonators. Fortunately, the addi- 
tional multiplier coefficients corresponding to this increased 
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bandwidth were insignificant i.e., = 0.000 and = 0.043. In 

fact, if is ignored and the resulting impulse response is compared 
with the ideal pulse, it has been observed that the maximum error 
remains at 0.84% 

3.3.3 Additional Comments 

Although the linear-programming technique for the design of 
digital filters has been stated as a time-domain procedure, the two 
examples made use of pulse bandwidth information for specification 
of R, the number of unknowns. The bandwidth, and consequently R, 
can always be obtained from a rigorous Fourier analysis of the- 
specified waveform or to a lesser degree of accuracy from an inspec- 
tion of the waveform shape. Alternatively, digital filters can be 
designed entirely in the time domain. For example, a matched filter 
for the Ricker seismic wavelet [68] has been designed without spec- 
tral information. The desired impulse response and the three 
attempts that were required to determine the final filter coeffi- 
cients are illustrated in Figs. 3^8 through 3^10. Constraint equa- 
tions and values used for the parameters N and T were defined from 
the Ricker wavelet plot to be 


h(.022) = -0.443 h(.032) = 0.198 

h ( . 028) =0 h’ (.032) < 0 


where it was desired to maximize 


(3.22) 


Z = h' (.032) 


(3.23) 




Filter Design (R-10). 
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with NT=0.44s and T = 0.00088s. The number of unknowns was not 
specified from the wavelet bandwidth but was determined by iteration 
If R was chosen too small, the linear-programming problem would not 
have a feasible solution, but if it was too large the resulting 
impulse response would have an unnecessarily high frequency content. 
After an initial estimate of ten for R, the impulse response of Fig. 
3-9 was obtained. Although the response satisfies the constraint 
equations, it is very evident from the high frequency envelope that 
R is too large. Consequently, R was reduced to seven in Fig. 
3-10(a), and finally to six, before the satisfactory impulse 
response of Fig. 3-10(b) was obtained. The corresponding changes 
in the set of multiplier coefficients {H^} can be observed from the 
results given with each figure. Objective function values were 
Z = 0 for each case. The application of linear programming to 
frequency-sampling digital filter design for baseband-channel 
equalization, which will be considered in the next chapter, employs 
only time-domain specifications. Therefore, an iterative procedure 
such as just discussed is used to determine a satisfactory solution. 

All examples and equations derived in this chapter have been 
for the design of frequency-sampling filters. No application of 
linear programming to pulse-shaping transversal filters has been 
considered, since the unknown multiplier coefficients for this fil- 
ter are specified simply by samples of the desired pulse shape taken 
at T-second intervals. Also, it is a very inefficient way to gen- 
erate pulses since a large number of multipliers are required, e.g., 
the raised cosine pulse of Fig. 3-4 would require 60 multipliers for 
a transversal filter as opposed to 16 for a frequency-sampling 
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filter. However, the design of transversal filters for channel 
equalization is neither a trivial problem nor an inefficient solu- 
tion, as will be illustrated in the next chapter. 



CHAPTER 4 


DIGITAL FILTER DESIGN FROM KNOWN CHANNEL SPECIFICATIONS 

The application of linear programming to the time-domain design 
of frequency-sampling digital filters has been discussed in Chapter 3. 
An extension of that technique will be developed in this chapter for 
the purpose of equalizing nonideal baseband data transmission channels 
to reduce intersymbol interference. Time-domain equations are 
derived for the received pulse at the channel output so that the pulse 
can be constrained to have a desired shape before sampling and thresh- 
old detection. Linear programming is used to design either a frequency- 
sampling or transversal digital filter which in cascade with the chan- 
nel will result in the desired pulse at the sampler when an impulse is 
applied to the equalizer. 

4.1 DIGITAL FILTER EQUALIZERS 

The technique of digital-filter equalization is illustrated in 
Fig. 4-1. The digital equalizer which replaces the normal trans- 
mitter filter will be designed to have a "predistorted impulse 
response". This predistortion is used to compensate for distortion 
that occurs during passage through the channel, thus resulting in an 
ideal pulse y(t) at the output. Each symbol a^ is considered to be 
the weight of an impulse that is applied to the digital filter to 
generate a predistorted pulse. The required D/A converter and any 
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receiver filter which would be used to reduce channel noise are 
both considered as part of the channel and must be included in 
c(t), the total channel impulse response. Equations for y(t) and 
y' (t) as functions of both filter multiplier coefficients and values 
of c(t) will be derived in this section for two finite impulse- 
response filters. The relationships between c(t) and c c (t) will 
also be defined. 


4.1.1 Frequency-Sampling Filter 

Equations (3.4) and (3.5) were derived in Chapter 3 for the 
impulse response of the frequency-sampling filter. By assuming that 
H^j /2 = 0 if N is even and redefining the upper summation limit to be 
R-l, the single equation 


h = h (mT) 
m 



I 


(-D 


2 \ 


N 


cos (2irkm/N) 


(4.1) 


is obtained for the predistorted pulse that is to be applied to the 
channel input. Since the sampling process is modeled with an impulse 
modulator in z-transform theory, the digital filter unit-impulse 
response can be rewritten as 

N-l 

h*(t) = h(mT) 8 (t - mT) . (4.2) 

m=0 

Consequently, the channel output for a. = 1 is given by 


y(t) = h*(t) * c(t) 


N-l 

= h(mT) 

m=0 


t 

J 6(t - mT) c(t - x)dx . 
0 


(4.3) 
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The limits of integration on the convolution integral can be verified 
with the graphical-convolution sketch of Fig. 4-2 (a). Applying the 
impulse-function sifting property to (4.3) yields 


y(t) = ^ h(mT) c(t - mT) 

m=0 


(4.4) 


where 


L = < 


I [ t/T ] 


N-l 


t < NT 


t > NT . 


(4.5) 


The channel output is obtained by substituting (4.1) into (4.4) i.e., 


H, 


L R-l 


k 2H k 


y (t) = ^ ^ c(t - mT) + ^ ^ (-1)^ cos (2irkm/N) c(t - mT) . 


m=0 


m=0 k=l 


(4.6) 


Hence, the output pulse slope is 




L R-l 


y ' (t) = ^ c' (t - mT) + ^ ^ (-1)^ — ^ cos (2iTkm/N) . c' (t - mT) 

(4.7) 


m=0 


m=0 k=l 


These last two equations are used to calculate constraint equation 
coefficients during linear-programming design of frequency-sampling 
equalizers. 

4.1.2 Transversal Filter 

It was shown previously by substituting (2.3) into (2.1) that 
the transversal-filter impulse response values are equal to the filter 
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multiplier coefficients. Consequently, the filter unit-impulse 
response can be written as 

N-l 

h*(t) = ) H 6(t - mT) , (4.8) 

m=0 

where the multiplier coefficients for the transversal filter are now 

denoted by H . Hence, the channel output is 
m 

L 

y (t) = h*(t) * c(t) = y (-l) m H c(t - mT) , (4.9) 

m=0 

where the parameter L was defined previously in (4.5). The factor 

(-l) m has been included because the set {H } will be the solution of 

m 

a linear-programming problem, and thus must be nonnegative. Without 
this factor, only nonnegative impulse response transversal filters 
would be designed, and for lowpass channels the output can not be 
constrained to have the desired zero crossings. Restricting every 
other multiplier coefficient to be negative simplifies the linear- 
programming problem of finding a feasible solution. The derivative 
of the output pulse is given by 

L 

y ' (t) =' ^ (-I)™ H m c'(t - mT) . (4.10) 

m=0 

Equations (4.9) and (4.10) are the desired expressions to be used 
during linear-programming designs of transversal equalizers. 
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4.1.3 Channel Impulse Response 

Since the derivative of the total channel impulse response 
appears in both (4.7) and (4.10), continuity requirements on c(t) 
must be investigated. The D/A converter will be modeled mathe- 
matically with a zero-order hold impulse response, i.e., 

g h (t) = u(t) - u(t - T) (4.11) 


where u(t) represents the unit-step function. Hence, if the receiver 
filter is considered to be included in the channel impulse response 
c c (t), the total channel response is given by 

c ( t) = c c (t) * g h (t) . (4.12) 


With the graphical convolution of Fig. 4-2(b), it can then be seen 


that 


c(t) 



c (x)dx 
c 


c (x)dx 
c 


0 < t < T 


t > T . 


(4.13) 


Therefore, it follows from Leibnitz's rule that the derivative exists 
and can be written as 


c'(t) =< 


c (t) 0 < t < T 

c — 


c (t) - c (t - T) t > T , 

c c 


(4.14) 


if c c (t) is continuous over the interval 0 < t < °°. The t = 0 point 
has been excluded in (4.13) and (4.14) because some channels of inter- 
est will have a step discontinuity at the origin. However, all channel 
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responses to be considered will be continuous for t > 0 and therefore 
(4.14) is valid. 

4.2 COMPARISON OF EQUALIZER DESIGNS 

Results of various design examples are presented in the remainder 
of this chapter to illustrate some of the distinctive characteristics 
of digital-filter channel equalization and guidelines for the use of 
linear programming to design these equalizers. The ease of design and 
relative performance of the frequency-sampling and transversal equal- 
izers are compared by using both techniques to equalize the channel 

c (t) = 4iT^t e ^ 7rt u ( t ) «-»- C (f) = » (4.15) 

c c (i + jf r 

Typically, the channel unit-impulse response will not be symmetric 
about the response peak. Consequently, unlike the frequency-sampling 
filter design examples of Section 3.3, both sides of the pulse must 
be constrained for equalizer designs. 

4.2.1 Frequency-Sampling Equalizer Design 

Since the channel described in (4.15) has a normalized break- 
point frequency f^ = 1 Hz, which is assumed to equal the Nyquist 
bandwidth, l/2t^, it follows that the minimum bit interval, i.e., 
the Nyquist interval, is t^ = 0.5 s. Consequently, the waveform 
specified for the equalized channel output was the raised-cosine 
pulse of Fig. 3.4 with t^ = 0.5 s, which resulted in the following 


set of constraints: 
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y(0. 5) = 0 y (1. 5) = 1 y(2.5) = 0 y'(l.O) ^0 

y(0. 75) = 0 y(2.0) = 0 y(3.0) = 0 y’(2.0) £0 (4.16) 

y(l.O) = 0 y(2.25) = 0 y’ (0.5) >_ 0 y' (2.5) <_ 0 . 

The objective function to be maximized is 


Z = -y ' (0.5) - y’ (1.0) + y’(2.0) + y'(2.5) , (4.17) 

which attempts to flatten the response on both sides of the pulse peak. 

The frequency-sampling equalizer design is also not as straight- 
forward as the pulse-shaping filter design in the determination of the 
parameters N, T, and R. The point NT/2 corresponding to the equalizer 
impulse response peak no longer coincides with the peak of the response 
being constrained, namely the channel output, since there is a delay 
in passage through the channel. Consequently, the time delay must be 
either estimated and subtracted from the time of the constrained peak 
to give the proper NT/2 value or it can be obtained by iteration. 

The delay was estimated from a sketch of the phase characteristic for 
this channel to be approximately 0.25 s resulting in the value NT = 

2.5 s. Thus, for an assumed sampling interval of T = 0.2^ cor “ 
responding number of impulse response samples is N = 25. The para- 
meter R is left as a variable to be determined by iteration. 

The impulse-response output of a known transmission channel, 
(4.15), which is in cascade with a frequency-sampling equalizer has 
now been constrained by (4.16) and (4.17) to have a desired impulse 
response. The equalizer consists of a comb filter with 25 delay 
elements and a set of R resonators with unknown multiplier coeffi- 


cients {H^}. The LINP0 subroutine can now be used to solve this 
linear-programming design problem to give the {H^}. However, the 
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coefficients of the unknown (H^.) as defined by (4.6) and (4.7) must be 
first calculated. Consequently, an auxiliary subroutine called C0CAL2, 
which is described in Appendix C, is used to prepare equalizer design 
problems for a LINPO solution, similar to the way C0CAL1 was used for 
pulse-shaping filter designs. 

The use of C0CAL2 and LINPO to solve (4.16) and (4.17) for the 
unknown is presented as an example in Section C.2. Results of that 
solution are illustrated in Fig. 4-3 where the 25 discrete values of 
the impulse response corresponding to the designed equalizer are 
plotted. Also shown is the output of the data channel when this 
impulse response is applied to a D/A converter at the channel input. 

By comparing the channel output with the constraint set of (4.16) it 
can be seen that all constraints are satisfied. Nevertheless, the 
pulse shape does not approach the desired characteristic since there 
is a large negative overshoot. Also, the channel output peak doesn't 
occur exactly at t = 1.5 s since the maximum channel delay is actually 
0.159 s rather than the 0.25 s used to calculate the product NT. Con- 
sequently, the equalizer parameters N, T, and R were varied in an 
attempt to improve the design but no significant improvement could be 
obtained over the illustrated design of N = 25, T = 0.1 s, and R = 12. 

The negative overshoot, which is a characteristic of frequency- 
sampling equalization, can be rationalized by sketching the graphical 
convolution of the equalizer output given in Fig. 4-3(a) with the 
channel impulse response as in Fig. 4.2(a). Considering just the six 
most significant values of h*(t), it can be seen that the negative 
overshoot is formed as c(t) translates past the two negative h*(t) 
values to the left of the h*(t) peak. If these two values were not 
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present, the y(t) pulse peak would be monotonically formed without an 
overshoot as c(t) translates past the two positive h*(t) values. Next, 
as desired the y(t) pulse would be immediately forced toward zero by 
the two negative values on the right of the h*(t) peak. Unfortu- 
nately, this desired nonsymmetrical impulse response cannot be 
obtained from a frequency-sampling filter with real multiplier coef- 
ficients, and thus the negative overshoot will always be present. 

4.2.2 Transversal Equalizer Design. 

Most of the design problems encountered during the frequency- 
sampling equalizer design are eliminated when transversal equalizers 
are considered. Since the transversal filter has the same number of 
impulse response values as multiplier coefficients, the parameters N 
and R are equal; thus, the design problem is simplified by reducing 
the number of unknowns that must be defined. Additionally, unlike 
the frequency-sampling equalizer, the LINPO requirement of real 
multiplier coefficients does not imply that the impulse response for 
this filter need be symmetrical. Since this filter also does not 
have a characteristic impulse-response peak at NT/2, there is no 
need to accurately estimate the delay through the channel to define 
the NT product. The discrete input to the channel is the set of 
multiplier coefficient values which are applied to the channel in 
sequence, one every T seconds. Consequently, the primary require- 
ment on the NT product is that the channel's response to the impulse 
response values in this time interval be of sufficient length to 
span the constrained output interval. 
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As an example, a transversal filter was designed to equalize the 
channel described by (4.15) so that the constraint set of (4.16) and 
(4.17) would be satisfied. All that remains to be specified for a 
LINPO solution are the values for N and T. The same parameters are 
used as in the previous example, namely a sampling interval T = 

0.2t^, and N = 25 impulse response values. Hence, the resulting NT 
product of 2.5 s is approximately equal to the constrained time 
interval and should be more than sufficient for the channel output 
to satisfy the constraints. It will be shown in Section 4.3 that 
the sampling interval has a direct effect on the transmitted energy 
requirement and can be optimized to reduce this energy. 

The C0CAL2 subroutine was used to calculate coefficients defined 
by (4.9) and (4.10), and then LINPO was employed to obtain a solution 
to the design problem. The solution, which is described in more 
detail in Section C.2, is illustrated in Fig. 4-4. The transversal 
equalizer requires only three multiplier coefficients as indicated 
by the three nonsymmetrical impulse response values in Fig. 4-4 (a). 
Since the last response value is at t = 1.8 s, the parameter N could 
have been reduced to 19 with no effect on the results. However, N 
will normally be set to a large value to eliminate the problem of 
choosing exactly its minimum possible value. By using this approach, 
the LINPO solution will indicate the true number of impulse response 
values required and specify the excess values to be zero. Thus, the 
multipliers and shift-register delays corresponding to these zero 
values need not be implemented, however, the channel response will 
occur earlier. The channel output produced by applying the three- 
valued impulse response to the channel is plotted in Fig. 4-4(b). 
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All constraint equations are satisfied and there are no overshoots to 
cause jitter problems. 


4.2.3 Timing Jitter Sensitivity 

The sensitivities to timing - jitter of the two types of digital 
equalizers were compared by considering the equalization of channels 
with more severe characteristics. Both frequency-sampling and trans- 
versal filters were designed to equalize the channels 


c 

c 


(t) 


( 2-rr) n n-1 -2fft 

(n - 1) ' 6 


u(t) 


n = 3, 4, 5 


(4.18) 


in the manner of the previous two design examples with the same con- 
straint set, T, N, and R values. Channel outputs, y(t), were then 
plotted so that the intersymbol interference (ISI) distortion could 
be calculated for each design as a function of timing jitter. This 
distortion is defined by 


D = £ |y(£t b + t p )| (4.19) 

£=0 

£^3 

where the time t is the percentage of the Nyquist interval that the 
sampling times at the receiver are displaced by jitter. The distor- 
tion was next expressed as a percentage of the pulse peak, y(3t, ) 

b 

for each design and plotted in Fig. 4-5 as a function of t . It is 
obvious that there is an appreciable advantage in using transversal- 
filter equalization for channels where jitter may be present. 
Furthermore, even if there was perfect timing, a frequency sampling 
equalizer could not be designed with LINPO for the n = 5 channel in 



ISI Distortion as Percent of Pulse Peak 
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Jitter as Percent of Nyqulst Interval 


Fig. 4-5. Comparison of Timing-Jitter Degradation. 
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(4.18). Since the transversal equalizer has the inherent advantages 
of being less sensitive to jitter and easier to design, the remainder 
of this chapter will be devoted to its design. 

4.3 OTHER TRANSVERSAL-EQUALIZER CONSIDERATIONS 

It has been shown that the only unknowns to be defined before 
using LINPO for a transversal-equalizer design are a set of constraint 
equations; N, the total number of impulse response values; and T, the 
sampling interval. The selection of these parameters will now be con- 
sidered. 

4.3.1 Effect of Constraint Set 

Only one constraint-equation set, (4.16), was used in the pre- 
vious section to design the transversal equalizer. Consequently, 
the effects of other constraint sets on the equalizer design were 
determined for the channel of (4.18) with n = 6, for which the 
distortion effects were more pronounced. Summarized in Table 4.1 
are nine of the sets of equations that were used. Starting with 
the set of seven constraints (NC =7) 

y(1.0) = 0 y(l. 5) = 1 y(2.0) = 0 

y(1.25)>0.5 y(1.75)>_0.5 y'(1.0)>_0 (4.20) 

y' (2.0) < o , 

additional constraints were added for each set, usually another time 
point at which y(t) was to be forced to zero. The parameters T ■= 

0.1 s and N = 40 were used for all designs. For complete elimina- 
tion of ISI distortion, y(t) must be forced to zero at every time 
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y(t) = o 
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point which is an integer multiple of t^ from the pulse peak. If 
only a finite number of such points are so constrained, some residual 
distortion will exist. This is illustrated by the transversal filter 
curve (n = 5) in Fig. 4-5, where for zero jitter a 0.2% residual dis- 
tortion still exists. Hence, the number of time points at which y(t) 
was constrained to be zero was increased until the residual distor- 
tion became insignificant. Typical output pulses resulting from the 
designs are plotted in Fig. 4-6, with other results summarized in 
Table 4.2. The tabulated results for NC = 11 through 15 indicate 
that an additional multiplier is required in the equalizer for each 


TABLE 4.2 


COMPARISON OF EQUALIZER CHARACTERISTICS 


Number of 
Constraints 

Multipliers 

Required 

Residual ISI 
Distortion 

Energy 

Transmitted 

7 

7 

1.240 

75.6 

11 

11 

1.739 

120.8 

12 

12 

0.222 

149.7 

13 

13 

0.083 

150.0 

14 

14 

0.035 

150.1 

15 

15 

0.015 

150.1 

16 

15 

0.015 

344.8 

17 

17 

0.002 

554.3 

21 

20 

0.000 

556.0 


additional time point constrained to zero at an integer multiple of 
t^ from the peak. However, the energy transmitted to the baseband 
channel, i.e., 



/ 


g 2 (t)dt = T 


N-l 


I 

m=0 



(4.21) 




Fig. 4-6. Pulse Shape Sensitivity to Constraints. 
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where g(t) = h*(t) * g^(t) approaches a constant value. The resid- 

ual distortion decreases with each added constraint and is approaching 
zero; consequently, the multiplier value h^ required to satisfy each 
added constraint also is approaching zero. Therefore, the summation 
of (A. 21) approaches a constant value. It can be seen from Table 4.1 
that the last three constraint sets, NC = 16 through 21, force zeros 
at time points between those resulting from the first six sets. These 
constraints result in the ISI distortion being reduced to an insig- 
nificant value at the cost of additional multipliers and a greatly 
increased transmitted energy. Attempting to force y(t) to be zero 
at time points this close together requires components of g(t) to be 
transmitted at higher frequencies where the channel attenuation is 
much greater. 

The binary data sequence 101011110010 was transmitted through 
the unequalized channel and then through the same channel using the 
NC = 15 constraint equalizer to illustrate the improvement obtained 
by using this design technique. Since the time between bits of the 
data sequence is 0.5 s and the sampling interval of the equalizer 
is 0.1 s, four equally-spaced zeros must be applied to the equalizer 
between each bit of the sequence so that the proper waveform will be 
generated. Comparison of the received waveforms in Fig. 4-7 shows 
that the equalized channel output has no intersymbol interference. 
Conversely, for the unequalized channel, detection errors are pos- 
sible even in the absence of noise. 

4.3.2 Effect of Sampling Interval 

It has been indicated previously that N, the number of impulse 
response values, should be made quite large so as to insure that 
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(a) No Equalization 
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Fig. 4-7. Effect of Equalization on Data 
Sequence 101011110010. 
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the product NT spans the constrained output time interval. Only the 
pertinent number of shift register stages will be retained in the 
final design. This leaves only the effect of the parameter T to be 
determined. 

The value of T is bounded on the lower end by the available 
logic speed and on the upper end by the Nyquist interval, i.e., 
t c _< T t^ , where t is the total time required for any one of the 
N multiplications plus one multi-input addition needed to calculate 
each transversal filter output sample. Hence, the effect of varying 
the parameter T between these limits was observed for equalizers 
designed with the NC = 15 constraint set and N = 40. The most sig- 
nificant result from these designs was the change in transmitted 
energy defined by (4. 21). These changes are displayed in Fig. 

4.8 for a number of channels. As T approaches t, in value there is 

b 

a great increase in required energy. A small percentage of the 
increase is a consequence of multiplying by T in (4.21), but the 
primary increase is caused by the reduction in bandwidth of the 
g(t) pulse applied to the channel. As T increases, the high 
frequency components of g(t) are severely attenuated as a result 
of the sin(irfT) /irfT frequency characteristic of the D/A converter. 
However, higher frequency components of g(t) must still be trans- 
mitted to satisfy the constraints on y(t). Hence, a very large 
increase in the transmitted energy. Conversely, as T approaches 
zero the bandwidth of the transmitted pulse becomes excessively 
large and since there is no energy minimization in the design 
constraints, these higher frequency components are not forced to 
zero, resulting in a slight increase in the transmitted energy. 
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It can be observed from the curves of Fig. 4-8 that a ratio of 

T/t^ ~ 0.4 minimizes the transmitted energy. For example, if the 

equalizer with NC = 15 and T = 0.1 s were redesigned using T = 0.2 s 

2 

the transmitted energy would be reduced from 150.1 V s as indicated 

2 

in Table 4.2 to approximately 100 Vs. It was also observed from 
the series of equalizer designs that the tails of the received 
pulses increased in magnitude as T varied from the minimum energy 
value of 0.4t^. The number of multiplier coefficients required did 
not change appreciably as T varied. 

4.3.3 Location of Equalizer 

It has been assumed throughout this chapter that the equalizer 
was positioned at the channel input, because there is a speed advan- 
tage in using this location for high data-rate applications. If it 
is assumed that 0.4t^ = T = t c then the maximum data rate through 
the equalized channel is 1/t^ = 0.4/t^. For the normal transversal 
filter with N multipliers, t c consists of the time required to 
perform one multiplication and one multi-input addition. However, 
with the equalizer at the channel input, multiplication is not 
required. Only one of two possible bits, a ONE or ZERO, will be 
applied to the equalizer input, and eventually to the equalizer 
multipliers, during each sampling interval. Therefore, each multi- 
plier can be replaced with a storage register which contains the 
multiplier coefficient; hence, t^ is the time required to read and 
add the appropriate coefficients. By reducing t c> the maximum data 
rate has been increased. Multipliers would normally be required 
for equalizers at the channel output, since there is a wide range 



93 


of possible inputs from the analog channel. However, there are 
advantages to equalization at the channel output. For example, in 
the case of the adaptive equalizer there would be no need to trans- 
mit the impulse responce c(t) back to the channel input in order to 
update the equalizer coefficients. If equalization is desired at 
the channel output, the linear programming design procedure of this 
chapter -is still valid. This is illustrated by comparing the chan- 
nel output of Fig. 4-6 (c) with the equalizer output plotted in Fig. 
4-9 (b). The former response was obtained with a transversal equal- 
izer designed and used at the channel input, while the latter is 
the output with the equalizer at the receiver. 







CHAPTER 5 


CONCLUSIONS AND RECOMMENDATIONS 

It has been shown that linear programming optimization tech- 
niques can be employed effectively in the time domain to design 
finite-duration impulse-response digital filters. It has been 
further demonstrated that this technique has direct application to 
the design of frequency-sampling filters to be used for pulse gen- 
eration or matched-filter detection, and both frequency-sampling 
and transversal filters to be used for equalization of transmission 
channels with known impulse responses. 

5.1 CONCLUSIONS 

Frequency-sampling filters, which have not been designed pre- 
viously in the time domain by any technique, can be defined readily 
by using linear programming to specify filter multiplier coefficient 
values. It was illustrated with examples that a simple design pro- 
cedure which can be employed is to: (1) constrain the impulse 

response equation of the filter to have desired values at critical 
time points, (2) define the number of impulse-response values, maxi- 
mum number of resonators and sampling-interval parameters for the 
filter, and (3) solve the resulting set of constraint equations with 
a linear programming algorithm. It was also shown that the number 
of response values and sampling-interval parameters are defined by 
the significant duration and point of symmetry of the filter impulse 
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response, whereas the number of resonators are defined by the impulse 
response bandwidth, if available, or by iteration. Comparisons of 
filter designs with desired pulses indicate that this technique yields 
very accurate pulse shapes with minimum hardware requirements. 

This design procedure can be extended readily to the design of 
digital equalizers for known transmission channels by constraining 
the output of an equalizer-channel cascade rather than a filter 
impulse response. It has been concluded by studying the equalization 
of typical transmission channels that the transversal equalizer is 
easier to design, requires considerable less hardware, and is less 
sensitive to timing jitter than the frequency-sampling equalizer. 

It was further observed from a study of transversal equalizer con- 
straint sets that an additional multiplier is required for each 
zero specified in the output response when spaced at Nyquist time 
intervals, but that the required transmission energy approaches a 
constant value. However, if zeros are forced in the response at 
closer time intervals, the transmitted energy requirements greatly 
increase. Results of varying the sampling interval in a series of 
equalizer designs indicate that this energy can be minimized by 
proper choice of the ratio of sampling to Nyquist intervals. 

The ability developed in this study to work directly from chan- 
nel time-domain characteristics allows the effect of both amplitude 
and phase distortions to be considered simultaneously. Unlike 
frequency-domain designs, which usually require separate filters to 
compensate for each source of distortion, one equalizer can be 
designed easily to compensate for both impairments. Furthermore, 
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the design technique is valid for equalizers positioned at either the 
transmitter or receiver. 

5.2 RECOMMENDATIONS FOR FURTHER STUDY 


Several possibilities exist for improvement in the equalizer 
design technique of this study. The feasibility of modifying or 
replacing the linear programming optimization algorithm should be 
investigated. For example, if an objective function in the form of 


Z = 


C 1 X 1 


C 2 X 2 


+ 


C R X R 


(5.1) 


was used, the transmitted energy (4.21) could be minimized directly. 
This would require a nonlinear programming routine because of the 
nonlinear objective function. Another possibility is to use integer 
linear programming, i.e. , the variables x^ of the problem are inte- 
gers. The integer variables would more accurately represent the 
possible multiplier coefficients, because the coefficients must be 
implemented with a finite number of bits. Thus, only a discrete, 
equally-spaced set of values would be possible solutions to the 
linear programming problem. 

The linear programming technique could be applied to the design 
of other classes of recursive digital filters. Since these filters 
typically have an infinite impulse response, it may be possible to 
minimize the residual intersymbol interference with fewer multiplier 
coefficients using less transmitted energy than the transversal 
equalizer requires. 

Finally, the use of the time-domain linear-programming technique 
to design filters should be investigated for other applications. For 
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example, the design of resonators for digital formant synthesizers 
would probably be a good application. Since the speech process can 
be described as the excitation of resonators by a series of impulses, 
specifications for synthesizer designs are often in the form of 
desired impulse responses. 


APPENDIX A 


LINEAR PROGRAMMING FUNDAMENTALS 

The general linear-programming problem is to solve for the R 
non-negative variables x. which maximize the objective function 


Z = 


c l*l + C 2 X 2 + 


+ C R X R > 


(A.l) 


subject to the M linear constraint equations 


< 

*il x l + a i2 X 2 + + a iR X R |>| b i 1 = 1,2 M 


(A. 2) 


One solution to this problem, namely the simplex procedure outlined 
in Section 3.2, is the subject of this appendix. Theoretical aspects 
of the simplex method are discussed first, followed by an example to 
illustrate the method. 

The simplex procedure involves a series of iterations, each of 
which generates a new basic feasible solution to (A. 2) that yields a 
greater value of the objective function than the previous solution. 
The procedure continues until a finite maximum value of the objective 
function is obtained or an unbounded solution is indicated. The 
aforementioned terms, e.g., basic feasible solution, were defined in 
Section 3.2. 

Any inequalities in (A. 2) are converted to equalities by the 
inclusion of non-negative slack and surplus variables. For example, 
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assuming in (A. 2) that holds for U equations, for V-U 

equations , and "=" for M-V equations then 


and 


R 


I Vj 

3-1 


+ X R+h ' b h 


R 


I Vi. - 

3-1 


X R+k ' b k 


R 


V a .X. = 

Z_< PI 3 


j-1 


h = 1 , 2 , . . . ,U , 


k = U+l,. . . ,V , 


(A. 3a) 


(A. 3b) 


p = V+l , . . . , M . 


(A. 3c) 


Hence, there are U slack variables x^ + ^ and V-U surplus variables x r+R ' 
The addition of these variables will have no effect on the final solu- 
tion since they will be added to the objective function with a zero 
cost, i.e., a coefficient c^ = 0, j = R+l, R+2,...,N. 

The M equations (A. 3) in N unknowns, N = R+V, can be written in 
vector notation* as 


A x = x^a^ + . . . + x^a^ = b_ , (A. 4) 


with the objective function written as 

Z = C 1 X 1 + ‘ • + c n x n = - - (A. 5) 

It will be assumed that the rank of A is M and that a solution to (A. 4) 


*A general matrix, e.g. , A will be denoted by a capital letter, 
whereas a single-row or single-column matrix, called a vector, will be 
denoted by a lower-case letter, e.g., x* A column vector will be 
denoted by parentheses, e.g., 2 C = (x^ ,x^ , . . . ,x^) and a row vector by 
brackets, e.g., c_ = [c^ , , . . . , c N ] . 
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exists, i.e., N >_ M. If N = M, there is a unique solution to the prob- 
lem and there is no need for linear programming. On the other hand if 
N > M there are a number of solutions and linear programming is used to 
determine the optimal solution. If either of the two assumptions does 
not hold, the simplex method will indicate that either there is no 
solution (due to inconsistent constraints) or there is redundancy in 
the constraint equations.* 

A. 1 ALGEBRAIC FORMULATION OF SIMPLEX PROCEDURE 

The simplex solution of the general linear-programming problem is 
an algebraic technique, which is best described in linear algebra 
terminology. Any set of M linearly independent columns of A forms a 
so-called "basis" for the vector space (E ) and the matrix ]} formed 
from these M columns is called the "basis matrix". Each column of _B 
is identical to a column of A, but the order of the columns in _B may 
be different from A. Corresponding to each basis matrix there exists 
a "basic feasible solution", x , to (A. 4), i.e., 

— D 

x g = B 1 b = (x B1 ,. . . ,x BM ) . (A. 6) 

The elements of the column vector x are called "basic variables" and 

— “D 

represent values for M of the N unknown variables x^ of (A. 4). The 
variable x^ which corresponds to the basic variable is determined from 
the column order of A in the basis matrix B_, e.g., if the third column 


*For a further discussion of inconsistency and redundancy see 
Section 4.6 of Hadley [67]. 
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of A is the first column of B then x = x_, etc. The remaining 

I>1 J 

N-M variables will be assigned zero values. 

For each basic solution, (A. 5) can be written as 


Z " C B1 X B1 + 


+ C BM X BM -B -B 


(A. 7) 


The component of c that corresponds to c . is again determined by the 
— Bl 

column order of A in B, e.g., c - = c_ if the third column of A is 
the first column of 15. The other N-M coefficients of c_ have no effect 
on (A. 7) since they are multiplied by zero-valued variables. The prob 
lem of selecting the columns of A which constitute the initial basis 
matrix will be discussed in Section A. 3. 

Each iteration of the simplex procedure consists of replacing one 
column of the basis matrix B. with one of the N-M columns of A not in 
the basis at that time. This substitution must result in a feasible 
basic solution (not just a basic solution which could contain negative 
variables) and increase the value of Z in (A. 7). Consequently, a pro- 
cedure will be discussed in the remainder of this section for deter- 
mining the column to remove from 15 so that a feasible solution is 
obtained and the column from A to substitute into _B so that Z is 
increased. However, certain variables must be defined first. 

Any column of A can be written as a linear combination of the 

til 

columns of 15 since B is a basis matrix. For example, the j column 
of A is 


a. = y 1 . b + ... + y„. = B y. . 

-J lj “I Mj -^1 — 


(A. 8) 
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The subscript order of the vector y__. components is important. The . 
first subscript refers to a column of B^ that the component is multi- 
plied by and the second refers to the column of A that is being 
expressed as a linear combination. Another function of the vector 
is to define a parameter, z_. , where 


z . = y . . c B1 + ... + y. . . c_._ = c D y. 
j ij B1 Mj BM — B 


(A.9) 


which is used to indicate possible increases in Z. 

It has been shown [67] that the column from A to be inserted 
into the basis should be a. where k is defined by* 


z, - c, = min {z. - c.; z. 
k he . 3 3 3 

and that the column of B to be removed is b 


■ c . < 0} 
3 

where . 


x. 


Br 

'rk 


= min 


Bi 

y ik 


y ik " ° 


> • 


(A. 10) 


(A. 11) 


These substitutions will insure a basic feasible solution after each 
iteration and will on the average give the largest possible increase 
in Z per iteration. The cost parameters, c^ in (A. 10), are constants 
fixed by the objective function, but the parameters z ^ must be recal- 
culated after each iteration. Likewise, the parameters x^ and y ^ 
of (A. 11) must be recalculated after each iteration. The repetitive 
calculations that are required for each iteration can be summarized 
as follows : 


*The set notation {aj; f(aj)} is used to define the set of 
elements aj which satisfy some restriction f(aj). Consequently, 
the element or elements in that set which have the smallest value 
are represented by min (a-i, f(a-O). 

3 j j 
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(1) Calculate the basic solution x with (A. 6) and the new basis 

— D 

matrix 15. 

(2) Calculate the vector £ with (A. 8) for every a^. not in the 
basis. 

(3) Calculate z^ with (A. 9) for every from step (2). 

(4) Calculate - c^ for every result of steps (2) and (3). 

(5) Select the column a^ to enter the basis with (A. 10). 

(6) Calculate the ratio x^/y^ for every y ^ > 0 of 

(7) Select the column b to leave the basis with (A. 11). 

— r 

(8) Calculate the new basis matrix by substituting a^ for b^ 

This procedure can end in one of two ways. If z. - c^ >_ 0 for all 
a. not in the basis, then the last basic feasible solution is an 
optimal solution. If there exists some column a.^. not in the basis 
such that z. - c, < 0 and y.. < 0 for all i = 1,2,...,M then an 

1 j ij “ 

unbounded solution exists. 


A. 2 EXAMPLE USING SIMPLEX PROCEDURE 

The example that was solved by graphical techniques in Section 
3.2 will now be worked by the simplex procedure to illustrate the 
steps of each iteration. The problem was to maximize 

Z = x i + x 2 

subject to the constraints 

x^ + x 2 2 , 

3x 1 + x 2 _< 15 , 
x^ + 2 x 2 _< 9 


and 
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After adding one surplus and two slack variables, the appropriate 
matrices are given by 

1 1-10 0 
A = 3 1 0 10 

1 2 0 0 1 

b = (2, 15, 9) , 

and 

c - [1, 1, 0, 0, 0] . 

As will be illustrated in Section A. 3 a basis matrix which will give 
an initial feasible basic solution is 



where the corresponding solution is 

x B = B _1 b = (2, 9,7). 

The vectors not in the basis are a_, ^ and a_^. Hence, from (A. 8) 

z 2 = B_1 1 2 = (1 » ” 2 ’ 1} 

and 

X 3 = B 1 a 3 = (-1, 3, 1) . 

) 

The solution is not unbounded since at least one > 0 for both 
vectors. From (A. 9), 

z 2 = £ B ^2 = [1 » °’ 0] = 1 

and 

z 3 ' % ^3 = _1 • 

Since 
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z 2 - c 2 = 1 - 1 = 0 , 

and 

z 3 - c 3 = -i - 0 = -i , 

it follows from (A. 10) that the vector to be entered is y^, i.e., 
k = 3. Since two £ components are positive, it follows from (A. 11) 
that two ratios must be calculated, i.e.. 



and 



which implies r = 2; hence b^ (° r equivalently a^) must be removed. 
Thus the second iteration starts with 


1 = = 


1-1 0 

3 0 0 

10 1 


Execution of the second repetition yields 


x B = (5, 3, 4) 

for the basic solution with 

2 2 = (1/3, - 2/3, 5/3) , 

£4 = (1/3, 1/3, - 1/3) , 
z 2 - c 2 = “2/3 , 

and 

z 4 “ c 4 » 

for the remaining variables. Consequently, the vector to be entered 
is given by k = 2, while the ratios 



107 


imply that b^ 


B1 

f 12 


x. 


= 15 


and 


B3 

32 


a_,_ is replaced with a^* 


12 

5 

The new basis matrix is 


1 = *y 2-2* 


and the corresponding basis solution is 


x = (4.2, 4.6, 2.4) . 
— B 


It then follows that 


and 


= (0.4, 0.2, -0.2), 

2 5 = (“0-2, 0.4, 0.6) , 

z. - c, = 0.2 , 

4 4 

z,. - c 5 = 0- 4. • 


Since z, - c. > 0 for the two vectors not in the basis, the basic 
1 1 ~ 

feasible solution is optimal. Hence, the complete solution is 


x = (4.2, 2.4, 4.6, 0, 0) , 

and the corresponding maximum value of the objective function is 

Z = jc x = 6 . 6 . 

A. 3 PROCEDURE FOR FINDING INITIAL BASIC SOLUTION . 

The initial basic feasible solution is obtained by always starting 
the procedure with an identity basis matrix, _I. If an M-by-M identity 
matrix can be obtained by manipulation of the columns in A, as would be 
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the case if all M constraint equations require slack variables, these 
M columns are used for the first basis matrix. Then 

Xg = B -1 b = I_ b = b . 

Since b^ is required to be non-negative, the initial basic solution 
is feasible. 

If a complete I matrix does not appear in A, a new set of con- 
straint equations are used so that an I matrix will be available, 
namely , 


A x + I_ x a = b^ . 


The vector 3 ^ represents artificial variables that are added to the 
original constraint equations with unity coefficients so that an 1_ 
matrix can be formed. The components of the artificial vector that 
correspond to constraint equations with slack variables will be zero. 

To illustrate the use of artificial variables consider starting 
the example in Section A. 2 without being given the initial basis 
matrix. The constraint equations with one artificial variable added 
are 


x x + x 2 - x 3 


+ x . = 2 
al 


and 


3x, + x„ + x, 

12 4 


x 1 + 2 x 2 + x c 


= 15 


= 9 


and the expanded coefficient matrix is 
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The initial basis matrix is 

®. — t ] ) 

where represents the artificial vector. Then the initial basic 
feasible solution is x = b. However, this is a solution to the new 

— D — 

constraint set, not the original set. Any solution to both sets must 
have 3£ a = 0. Hence, the iterations must drive the nonzero components 
of x^ to zero. This can be done with a two-phase objective equation. 

During the first phase, the artificial variables are assigned a 
cost of -1 and the other variables a cost of zero. Thus, the first 
phase objective function 


M 



i=l 

will be maximized by the simplex procedure. Since the simplex 
procedure prevents variables from becoming negative, the maximum 
value of Z* is zero and all artificial variables will be zero when 
this maximum is reached. At this point, which is the end of Phase 1, 
a basic feasible solution to the original problem is available. If 
the maximum value of Z* is negative, the artificial variables cannot 
be driven to zero and no feasible solution to the problem exists. 

The basic feasible solution which exists at the end of Phase 1, 
where max {Z*} = 0, is not optimal. Therefore, Phase 2 consists of 
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assigning the original costs to the legitimate variables and costs of 
zero to the artificial variables. The simplex procedure then obtains 
the optimal solutions in the normal manner. 

Continuing Phase 1 for the example, the columns of A not in the 
first basis are a,, a , and a. . Since c = (-1, 0, 0), it follows 
that 


= (1, 3, 1) , 

Z 2 = 0-. !. 2 > » 

and 

X 3 = (-1, 0, 0) . 

Then, since all non- artificial variables have zero cost values, 



Hence, k = 1 and the ratios are 


B1 

'll- 


= 2 , 


B2 

'21 


= 5, 


and 



Since these ratios imply that r = 1, column b_^ = e_^ is removed from 
the basis and replaced with _a . This gives 


B = 


10 0 
3 10 


and 


-B = (X 1 J X 4’ X 5 ) 


1 


0 


1 



Ill 


The artificial variable has been removed from the basis. Since any 
variable not in the basis is assigned a zero value it follows that 

3 

Z* = - ) x . = 0. 

Z_ i ai 

i=l 

Phase 1 has ended with a basic feasible solution. The resulting 15 
matrix is the one that was initially assumed in Section A. 2 to work 
the original example, i.e., Phase 2 of the simplex procedure. 
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SUBROUTINE LINPO 

LINPO is a double-precision computer subroutine written in Fortran 
IV which uses the simplex method of Appendix A to solve the general 
linear-programming problem. LINPO is a mnemonic for LINe ar JProgram- 
ming Optimization. This subroutine was developed by modifying an 
existing routine [69] so that the internal calculations of the iter- 
ative procedure could be followed and the routine could be used to 
design digital filters. Instructions for using the subroutine, an 
example solution, a flow-chart, and a program listing are contained 
in this appendix. 

B.l USER INSTRUCTIONS 


The call statement for the subroutine is CALL LINPO(A,B,IE,NEQ, 
NVA,NSP ,NT) with the required input parameters described in Table B.l. 


TABLE B.l 

LINPO INPUT PARAMETERS 

A: A two-dimensional double-precision array of real numbers repre- 

senting the objective-function and constraint-equation coeffi- 
cients. Coefficients of the objective function multiplied by 
-1 must be placed in the first row with each of the remaining 
rows containing the coefficients for one constraint equation. 

B: A double-precision array of non-negative numbers representing 

the right-hand sides of the constraint equations and entered 
in the same order as the A array. The first value, which cor- 
responds to the objective function, must be zero. 
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IE: An array of integer constants representing the type of inequality 

for each equation. The constant is one of the integers zero 
through three depending on the type, namely, objective function 
(IE = 0), less than or equal (IE =1), equal to (IE =2), or 
greater than or equal (IE = 3). The integers are entered in the 
same order as the B array. 

NEQ: Integer constant M + 1 equal to the number of constraint equa- 
tions including the objective function. 

NVA: Integer constant R equal to number of unknowns in the original 
constraint equations. 

NSP: Integer constant specifying tableau output, i.e. , tableau out- 
put (NSP =1), or no tableau output (NSP = 0). 

NT: Integer constant specifying type of digital filter to be 

designed, namely, frequency-sampling filter (NT = 0,1,2, 3) or 
transversal filter (NT = 4, 5, 6, 7). Furthermore, the LINPO sub- 
routine can be used for general linear-programming applications 
other than the design of digital filters by setting NT >_ 8. 

B. 2 EXAMPLE PROBLEM 


The problem chosen to illustrate the use of LINPO is the one 
which has been worked by a graphical technique in Section 3. 2 and 
by hand computation in Section A. 2, i.e., to find the maximum value 
of Z where 


Z = x^ + x 2 , 

subject to the constraints 


and 


x. 

+ 

x 0 > 2 

1 


2 — 

3x. 

+ 

x„ < 15 

1 


2 - • 

x i 

+ 

2x < 9 

1 


L — 


will now be solved using subroutine LINPO. The main program required 
is illustrated in Fig. B-l with the corresponding subroutine output 


given in Fig. B-2. 
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DOUBLE PRECISION A(27,50), B(4) 

DIMENSION IE (4) 

DATA A/-1. DO, 1. DO, 3. DO, 1. DO, 23*0. DO ,-1. DO ,1. DO, 1. DO, 2. DO, 1319*0./ 
DATA B/0. DO, 2. DO, 15. DO, 9. DO/ 

DATA IE/0,3,1,1/ 

DATA NEQ ,NVA,NSP ,NT/4 , 2 , 1 , 8/ 

CALL LINPO ( A , B , IE , NEQ , NVA , N SP , NT ) 

STOP 

END 


Fig. B-l. Main Program for LINPO Solution. 


Information is output in a tableau form, where the tableau 

entries are defined in Table B.2. Two rows of the tableau contain 

values for z. - c.. For Phase 1 of the simplex method, z. - c. 

J J J J 

values in the last row are used to determine the vector a^ to enter 
the basis, during which time the first row values are not significant. 
Conversely, for Phase 2 the first row values are used and the last row 
has no significance. The current objective function value is the last 
entry in the last row for Phase 1 and the last entry in the first row 
for Phase 2. The components of each ^ and the basic solution compo- 
nents X^, i = 1,2,...,M, are contained in the other M rows of the 


TABLE B. 2 
TABLEAU FORMAT 


A1 


A5 

XB 

z r c i 

• • • 

Z 5" C 5 

Z 

y n 


y 15 

X B1 

y 21 


y 25 

X B2 

y 31 


y 35 

X B3 

z r c i 

• • • 

Z 5" C 5 

z* 



VECTOR ENTERED A 0 ' VECTOR REMOVED 
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c o c c oc.cc. 

O O c c ccoc 

OCCC cccc 

ccoc cccc 

acccc cccccc 

XCCCC. * C C C C 

OCCC oooc 

CCCC' 

^ m ^ O ^ nC 


o o o o o 


© o o o o 




Fig. B-2. Tableau Output of Subroutine LINPO. 
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tableau. Also printed after each iteration are the vectors currently 
in the basis and the vectors which were entered and removed to form 
that basis. However, the initial vector entered and removed statements 
have no significance. All artificial vectors are represented by QO. 

Intermediate computer solutions contained in the tableau of Fig. 
B-2 agree with the results obtained by hand computation in Appendix A. 
for example the last iteration values in Section A. 2 were 

x B = (4.2, 4.6, 2.4) 

= (0.4, 0.2, -0.2) 

X 5 = <-0.2, 0.4, 0.6) , 

z 4 " c 4 = °.2 , 

z 5 - c 5 = 0.4 , 

and z - 6.6 

Each of these values can be found in the appropriate columns of 
Tableau 2 of Phase 2. 

B. 3 FLOWCHART AND PROGRAM LISTING 

A flowchart and program listing for the LINPO subroutine are 
given in Figs. B-3 and B-4, respectively. The calculations (A. 6) 
through (A. 11) , as discussed in Appendix A, make up the right-hand 
loop in the flow chart and each passage through the loop represents 
one iteration. The results of these calculations can be printed out 
at the end of each iteration in the tableau form. The left-hand 
loop determines if the solution at the end of Phase 1 is feasible 
and if so, sets up the problem for Phase 2. At the end of Phase 2 
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a branch of this loop calculates (from the optimal solution) the 

digital filter coefficient values Hq and (-1) 2H^, k = 1,2, . . . ,NVA-1, 

for a frequency-sampling filter design or (-l)™!! , m = 0 , 1, . . . ,NVA-1 , 

A m 

for a transversal filter design. These coefficient values along with 
the optimal solution are then printed. Filter coefficients for a 
frequency-sampling and a transversal filter design example are printed 
in Figs. C-6 and C-7 respectively. 



Function 
Set k^t 
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Pig. B-3. Flowchart for Subroutine LINPO 
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SUBROUTINE L IMP OH A , R, I E , K r: n , M VA , M SP , NT ) 

r. 

C THIS SUBROUTINE SOLVES THE GENERAL LINEAR-PROGRAMMING 

C PROBLEM BY THE SIMPLEX METHOD. THE SUBROUTINE. ACCEPTS 

C UP TO 25 CONSTRAINT EOUATIOMS AND AO VARIABLES. 

C 

DIMENSION L ( 27 ) , I E (NEO ) 

DOUBLE PRECISION A(27, 50), W( 2 7 ) , X , XMI N , B ( NED ) , H ( AO ) 
COMMON /HA' A ME / H 

INTEGER A0AI200), ADA 1 (200) ,A0, AN,BN, XN, BL, ON 
K J 11=0 

DATA BL/' •/ 

DATA XN/' X'/ 

DATA AN/ 'A'/ 

DATA BM/'R'/ 

DATA ON/'O'/ .• 

GO TO A22 

C WRITE TABLEAU HEADINGS AND VALUES (THRU ST ATFM ENT 42 I ) 

AOO TF (NSP.EO.O) GO TO A2 1 
JJJ1=JJJ 

TF(JJJ.GT.6) JJJ1-6 
K PH= 1 

TF(K.EO.l) K P H = 2 

IF(K.GT.l) KPH=1 

WR1TE( 6,401 )KKK ,KPH , 

401 FORMAT! 'O' ,51X, 'TABLEAU' ,12, IX, 'PHASE', 12 ) 

IF (KJ11.E0.0) GO TO A02 

"RI TE ( A, A03)L ( JK ) ,KJ1 1 
GO TO 405 

402 WRTTE(6,404) L(JK),KJ11 

403 FORMAT ( *0 ' ,20X, ' VECTOR ENTERED A ' , 1 2 ,4 IX , • V EC TOR. . 

I ' REMOVED A » , I 2 ) 

40^ FORMAT ( '0 ' ,20X, 'VECTOR EN T ERE D A ' , I 2 , 4 IX , • V EC TOR • , 

1 ' REMOVED O' , 12) 

405 WR I T E ( 6 ,406 ) 

40ft FORMAT! *0', ' BA S IS . VECTORS.' ) 

no 407 I=1,JJJ 
A OA 1 ( I ) =BL 
AOA ( I ) =AN 

I F ( I . EO.JJJ ) AOA 1 ( I ) =XN 

407 IF ( I .EO.JJJ ) AO A ( 1 ) =BN 
I F( JJJ.LE.6) GO TO 408 

WR ITE( 6,409 ) ( A 0 A 1 ( I ) , AOA ( I ) , I , I = 1 , J J J 1 ) 

GO TO 410 

408 JJJJJ=JJJ1-1 

W R I T E ( 6 , 40 9 ) ( AOA I ( I ) ,AOA( I) , I , 1 = 1 , JJJJ J ), 

1A0A1 ( J J J 1 ) ,AOA( JJJ1 ) 

40<? FORMAT! » + * ,24X,6( 2A1 , 12, 13X) ) 


Fig. B-4. Program Listing for Subroutine LINPO: 



410 CONTINUE 
JJ1=JJJ 

I F( JJJ.GT.6) J J 1 = 6 
DO 4111 = 1, III 
L LL = B L 
AO = BL 

IF(I.GT.l) LLL =L ( I ) 

IF ( I .GT.l ) AO=A N 
I F(L( I) .EO.O.AND. I.GT.l ) AO=ON 
IF ( I .EO. Ill) LL L =BL 
IF(I.EO.III) AO=BL 

411 WRITE! 6,412 ) AO , LL L , ( A ( I , J ) , J = 1 , J J 1 ) 

412 FORMAT! • ', 4X , A 1 , 12 , 1 1 X ,6 ( D 1 5 . 8 , 2 X ) ) 

K1K = 2 

J J 1 1 = 6 
JJ1 = JJ1 + 1 
J JJ 1 = JJ J-6 

413 IF ( JJJ 1 )421 ,421 ,414 

414 J J 1 1 = J J J 1+ J J 1 1 

IF ( JJJ1.GT.6) J J 1 1 =6 5 ! ; K IK 
I F( JJJ1 .LE.6) GO TO 415 

WR ITE! 6,417 ) ! A0A1 ( I ) ,AOA ( I ) , I , I = J J 1 , J J 1 1 ) 

GO TO 418 

415 JJ111=JJ11-1 

I F( JJ1.GT. JJ111) GO TO 416 

WRITE! 6,417) (A OAK 1 ) , AOA ( I ) , I , I = J J 1 , J J 1 1 1 ) , 

1 A OA 1 ( JJ 1 1 ) , AOA ( J J 1 1 ) 

GO TO 418 

416 WRITE(6,417) AOA 1 ( JJ 1 1 ) , AOA ( J J 1 1 ) 

417 FORMAT ( *0* ,23X,6(2A1, 12, 13X ) ) 

418 DO 419 1=1,111 

419 WR ITE! 6,420) ( A! I ,J ) ,J=JJ1 ,JJ11 ) 

420 FORMAT!' • , 1 8X , 6 ( D1 5 . 8 , 2 X ) ) 

JJ J1=J JJ1-6 

J J1=JJ1 1+1 
K1K=K1K+1 
GO TO 413 

421 CONTINUE 

I F ( K ) 457,456,457 

422 CONTINUE 

423 FORMAT! 47X, 14, 020,8) 

424 FORMAT! 'O' ,56X, 'FEASIBLE' ) 

425 FORMAT! 'O' ,56X, ' INFEASIBLE' ) 

426 F0RMAT(46X, 'VARIABLE' ,7X, 'VALUF' ) 

427 FORMAT! IX,//, 45X, 9H0PTIMUM Z,D20.8//) 

428 FORMAT (53X, 'BASIC SOLUTION'/) 

429 FORMAT! IX, 55 X, 'UNBOUNDED'// / ) 

430 FORMAT ( 1H , 'ROWS EXCEED ALLOWABLE MAXIMUM OF 25 


Fig. B-4 (Continued). Program Listing for Subroutine LINPO. 
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1* EQUATIONS'// 1 PROCESSING TERMINATED') 

431 F OR MA T ( I H , 'COLUMNS EXCEED ALLOWABLE MAXIMUM OF', 

1' 40 VARIABLES' //• PROCESS ING TERMINATED') 

C READ INPUT DATA (THRU STATEMENT 447) 

1 1 =N EQ 
J J = NVA+ 1 

432 IF ( I 1-25)434,434,433 

433 WRITE (6,430) 

GO TO 507 

434 I F( J J- 41)436,436,435 

435 WRITE( 6,431) 

GO TO 507 

436 II 1=11+1 
JJJ=JJ+II 

DO 437 1=1,111 
w ( I ) =0. 0 
L ( I ) =0 

DO 437 J=JJ,JJJ 

437 A ( I , J ) =0 . 0 

DO 438 J=1 ,NVA 
43 R A(III,J)=0.0 
JJJ=JJ 

DO 445 1=1,11 
J K= I E ( I ) 

439 IF ( JK ) 507 ,445,440 

440 IF(JK-3)441, 44 1,507 

441 JK =JK-2 
IF(JK)442, 444, 443 

442 A(I,JJJ)=1. 

L ( I ) = J J J 
JJJ=JJJ+1 
GO TO 445 

443 A( I , J J J ) = — 1 « 

JJJ=JJJ+1 

444 L ( I ) =0 

445 CONTINUE 

DO 446 1=1,11 

446 A ( I , J J J ) =B { I ) 

JJ=JJJ 

C WRITE DATA LOADED 

WRITE! 6,447 ) 

447 FORMAT( • 1 • ,54X, 'DATA LOADED') 

C START OF SOLUTION PHASE 

448 K KK=0 
JK =1 
K =0 

C ADD ARTIFICIAL VARIABLES TO PROBLEM (THRU STATEMENT 455) 

449 1=1 


Fig. B-4 (Continued). Program Listing for Subroutine LINPO. 
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450 1=1+1 

IF ( I-I II ) 45 1 ,455,455 

451 IF(L( I) 1450,452,450 

452 DO 454 J=1,JJ 

I F ( A ( I ,J ) >453,454,453 

453 A( I I I , J) = A( I I I , J )-A (I , J ) 

454 CONTINUE 
GO TO 450 

45 5 GO TO 400 

455 K= III 

C DETERMINE VECTOR TO ENTER BASIS (THRU STATEMENT 463) 

457 J= 0 

W ( K ) = 0 . 0 
L ( K ) =0 
45 B J=J + 1 

IF (J-JJ >459,462 ,462 

459 I F( DABS ( A(K, J ) ).L E.O. ID-08 ) A(K,J)=0. 

IF ( A (K ,J ) ) 460 ,45 B , 458 

460 TF(W(K)— A(K,J) )458,458,461 

461 W ( K ) = A ( K , J ) 

L (K ) = J 

GO TO 458 

462 I F ( L ( K ) ) 463,492,463 

463 K J =L ( K ) 

C DETERMINE IF SOLUTION IS UNBOUNDED 

DO 464 1=2,11 
I F ( A ( I , K J ) >464,464,465 

464 CONTINUE 

C WRITE UNBOUNDED 

WRITE( 6,429) 

GO TO 507 

C DETERMINE VECTOR TO LEAVE BASIS (THRU STATEMENT 473) 

465 1=1 
JK =0 

466 1=1+1 

IF ( I-I I )467 ,467 ,473 

467 I F ( A ( I, KJ) >466,466,468 

468 IF (DABS (A ( I ,KJ) ) -0 . 1 D-08 ) 469 , 469 , 470 

469 A ( I , K J ) =0. 

GO TO 466 

470 X=A( I ,JJ)/A( I , KJ ) 

IF (JK) 471 ,472,471 

471 I F(X-XMIN)472,466,466 

472 XM IN=X 
J K= I 

GO TO 466 

473 X =A ( J K , K J ) 

KJ 1 1 =L ( JK ) 


Fig. B-4 (Continued). Program Listing for Subroutine LINPO. 



L ! JK ) = K J 

C CALCULATE MEW TABLEAU VALUES (THRU STATEMENT 491) 

DO 474 1=1,111 

474 W (I ) = A ( I ,K J ) 

1 J = JK- 1 

DO 482 1=1 , I J 
DO 482 J=1,JJ 
I F ( A ( JK , J ) >475,482,475 

475 IF (DABS(A (JK,J) ) -0 . 1D-0 8 ) 47 6 , 476 , 477 

47 6 A ( JK , J ) =0. 

GO TO 482 

477 I F(W( I ) >478,482,478 

47 8 IF (DABS (W ( I ) ) -0 . ID-08 ) 479 , 479 , 480 

479 W ( i ) = 0 . 

GO TO 482 

480 A(I,J)=A(I,J)-W(I)*(A(JK,J)/X) 

IF ( DABS (A ( I , J ) >-0.10-0 8)481 ,481,482 

481 A ( I , J ) = 0 . 

482 CONTINUE 
I J= JK+1 

DO 490 I=IJ,III 

DO 490 J=1,JJ 

IF (A ( JK,J ) >483,490,483 

48 3 I F( DABS ( A( JK , J ) )-0.1D-08 >484,484,485 

484 A(JK,J)=0. 

GO TO 490 

485 IF(W(I ) >486,490,486 

486 IF(DABS(W( I) )-0.1 D-0 8 >487,487,488 

487 W(I)=0. 

GO TO 490 

48 8 A( I,J)=A( I,J)-W(I)*(A{ JK,J)/X) 

I F( DABS (A( I,J) >-0.10-08)489,489 ,490 

489 A ( I , J ) =0. 

490 CONTINUE 

DO 491 J=1,JJ 

491 A { JK, J) =A( JK, J )/X 
KKK = KK K + l 

GO TO 400 

492 IF (K-l >497,497,493 

C DETERMINE IF SOLUTION IS FEASIBLE 

493 IJ=JJ-1 

IF ( DAB S ( A ( K , J J ) >-0.10-0 8)494,4 94,496 

494 CONTINUE 

C WRITE FEASIBLF 

WRITE (6,424) 

DO 495 J = 1 , J J . 

495 A ( III, J)=0.0 
K = 1 


Fig. B-4 (Continued). Program Listing for Subroutine LINPO. 
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C 

496 
C 

497 


C 

498 


499 

500 
C 


501 

502 

503 


504 
C 

505 

506 

507 


K KK =0 
00 TO 457 
WRITE INFEASIBLE 
WRITE ( 6,425) 

00 TO 507 

WRITE OPTIMUM VALUE OF OBJECTIVE FUNCTION 
WRITE (6,427) All, JJ) 

WRITE (6,428) 

WRITE (6,426) 

WRITE LIST OF REAL VARIABLES 
00 498 1=2,11 

WRITE (6,423) L(I),A(I,JJ) 

I F (NT.GE.8) 00 TO 507 
WRITE! 6,499 ) 

FORMAT! Ill ,47X, 'DIGITAL FILTER COEFF I Cl ENTS 1 ) 

WRITE! 6,500) 

F0RMAT( '0* ,47X, 'STAGE* ,8X, 'VALUE' ) 

CALCULATE DIGITAL F ILTER ' COEF F IC I ENTS (THRU STATEMENT 504) 
DO 502 1=1 ,NVA 
DO 501 J = 2 , I I 

I F( L ( J) . EQ. I ) H( I ) = A (J,JJ) 

IF (L (J ).E0. I ) GO TO 502 
IF(J.EO.II) H(I)=0.000 
CONTINUE 
CONTINUE 

DO 503' 1=2, NVA, 2 
H ( I )= (-1 .DO) *( H( I ) ) 

IF (NT.GE.4) GO TO 505 
DO 504 1=2, NVA 
H ( I ) = ( 2 . 0 DO ) * ( H ( I ) ) 

WRITE LIST OF FILTER COEFFICIENTS 
DO 506 1=1, NVA 
WRITE (6, 42 3) I , H ( I ) 

RETURN 

END 


Fig. B-4 (Continued). Program Listing for Subroutine LINPO. 



APPENDIX C 


AUXILIARY SUBPROGRAMS 

The linear-programming subroutine LINPO used for the design of 
digital filters requires that coefficient arrays be calculated before 
the subroutine can be used. Since this procedure can become rather 
tedious due to the number and type of calculations required, especially 
for the channel equalization application, auxiliary Fortran IV sub- 
programs C0CAL1 and C0CAL2 were developed, where COCAL is a mneomonic 
for Coefficient CAL culator. C0CAL1 is a subroutine for the design of 
frequency-sampling filters which are to be used either for generating 
pulses to transmit data over ideal channels or for matched filter 
detection. C0CAL2 is a subroutine for the design of either trans- 
versal or frequency-sampling filters which are to be used to equalize 
nonideal channels. The nonideal channel must be defined by its unit- 
impulse response in a function subprogram. Consequently, a number of 
impulse-response equations corresponding to common data transmission 
channels have been programmed and are incorporated in a function sub- 
program called CC. Also included is a subroutine CHECK to be used 
in conjunction with C0CAL2 and CC. It calculates the unit-impulse 
response of the equalized channel after the design of the digital- 
filter equalizer. 
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C.l SUBROUTINE COCAL1 


The frequency-sampling filter impulse-response envelope was 
derived in Section 3.1 and expressed in (3.6) or (3.7) as 


h(t) 



E 


(-D 


2H k 


cos (27rkt/NT) , 


(C.l) 


where the upper limit on the sum has been replaced by R-l and the fre- 
quency sample has been assumed to be zero if N is even. The 

derivative of (C.l) is 


h'(t) 


4ir 

2 

NT 


R-l 

k(-l) k+1 sin (2-irkt/NT) . 

k=l 


(C. 2) 


Coefficients of the terms must be calculated for the objective 
function and every constraint equation of any LINPO frequency-sampling 
filter design. Also, for any LINPO design the constraint equations 
must be expressed in a form such that all b^ as defined in (3.11) are 
nonnegative. Furthermore, the objective function must be expressed 
as a maximization problem. C0CAL1 is a double-precision computer sub- 
routine that takes a more general linear-programming design problem 
and sets up the proper LINPO input arrays to satisfy these restric- 
tions. It also calls the LINPO subroutine to solve the problem. 

There are a variety of problems that C0CAL1 will accept. The 
objective function can be in any of the following forms: 


E a. h (t . ) , (C. 3a) 

J J 

j 


maximize Z = 
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or 


maximize Z = 


minimize Z = 


minimize Z = 


Y. “j h,( v • 

j 

T. “i h<t i> • 


I 


a. h' (t .) , 
3 3 


(C. 3b) 
(C. 3c) 
(C. 3d) 


where each a is either ±1. As many equations of one type as desired 

can be combined to form Z. Each individual constraint equation can 

be of the h(t) or h' (t) type with either positive, negative, or zero 

b. . 
i 


C.1.1 User Instructions 

The call statement for the subroutine is CALL C0CAL1 (CT,NTYPE,B , 
IE ,NOBJ ,N ,T ,NEQ ,NVA,NSP ,NT) where the input parameters IE, NEQ, NVA, 
and NSP were defined previously in Table B.l. The other input para- 
meters are defined in Table C.l. The parameter NVA was defined as 
the number of unknowns in the constraint equations and equaled R in 
the constraint equation formulation of Appendix A. For the low-pass 
and band-pass filter design problems of Chapters 3 and 4, R and con- 
sequently NVA, is equal to the number of unknown multiplier coeffi- 
cients to be determined optimally. The value for NT is generated 
internally in C0CAL1 and is an output parameter. 
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TABLE C.l 

C0CAL1 INPUT PARAMETERS 


CT: A double-precision array of real numbers representing critical 

time points. All time points at which the impulse response is 
to be constrained must be included. Furthermore, if more than 
one constraint equation is written at the same critical time, 
that time point should be repeated for each of these constraint 
equations. However, time points that appear in the objective 
function equation should not be repeated if constraint equa- 
tions have also been written at those points, as will normally 
be the case. If constraint equations have not been written at 
some of the objective function time points, those times should 
be included as the last members of the CT array. 

NTYPE: An array of integer constants representing the objective func- 

tion and constraint equation types. The first value in the 
array represents the type of objective function, i.e., "0" 
for an objective function in the form of (C.3a), "1" for 
(C.3b), "2" for (C.3c) and "3" for (C.3d). The remaining 
values in the array are in one-to-one correspondence with the 
constraint equation time points of CT and represent the type 
of each constraint equation, i.e., "0" for an h(t) equation 
and "1" for an h' (t) equation. 

B: A double-precision array of positive or negative real numbers 

representing the right-hand sides of the constraint equations 
and entered in the same order as the NTYPE array. The first 
value, which corresponds to the objective function, must be 
zero . 

NOBJ : An array of positive and negative integer constants defining 

the objective function. The first value is a positive inte- 
ger equal to the number of equations to be combined to form 
the objective function. The signs of the remaining integers 
define the a. of (C.3) while the associated absolute values 
define the locations in the CT array of the corresponding t ^ . 

N: A double-precision constant equal to the number of impulse 

response values. 

T: A double-precision constant equal to the sampling interval. 
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C.1.2 Example Problem 

The set of constraint equations that were defined in (3.16) for 
the raised-cosine design problem were 

h(3t, ) = 1 h(5t.) = 0 h'(4t,) < 0 

b b b — 

h(4t, ) = 0 ' h(6t, ) = 0 h'(5t.) < 0 

b b b — 

with the objective function 

Z = h' (4t, ) + h’ (5t, ) 
b b 


to be maximized. A main program is given in Fig. C-l as an example of 
the input required for a C0CAL1 solution of this problem. The corre- 
sponding output of the LINPO subroutine is illustrated in Fig. C-2. 


Since NSP = 0 for this example, none of the tableaus were printed 
during the linear-programming solution. The results printed first 


are the optimum value of Z and the optimal basic solution. The six 


values of the basic solution represent the unknown H^, k = 0,...,5. 


These values were used to calculate the digital filter coefficients 


H q and 2(1-1)^, k = 1,. 


,5 which are printed out last. 


DIMENSION NTYPE(7) , IE (7) ,NOBJ (3) 

REAL *8 CT(6) ,B(7) ,N,T 

DATA CT/3.0D0 , 4. 0D0,5 .0D0, 6 . 0D0 ,4. 0D0,5 . 0D0/ 

DATA NTYPE/1 ,0,0, 0,0, 1,1/ 

DATA B/0. 0D0 , 1. ODO ,5*0. ODO/ 

DATA IE/0,2,2,2,2,1,1,/ 

DATA NOBJ/2 ,2,3/ 

DATA N,T,NEQ,NVA,NSP/ 60. ODO ,0. 1D0 ,7,6,0/ 

CALL CO CALI (CT,NTYPE , B , IE.NOBJ ,N ,T ,NEQ ,NVA,NSP ,NT) 

STOP 

END 


Fig. C-l. Main Program for C0CAL1 Solution. 
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DATA LOADED 
FEASI BLE 


OPTIMUM 7. 0.0 



BASIC SOLUTION 


VARIABLE 

VALUE 


5 


0. 3333333 6D 

01 

1 


0. 10000000D 

02 

4 


0 . 50000000D 

01 

3 


0. 66666664D 

01 

6 


0.1 666667 3D 

01 

2 


0. 83333327D 

01 


DIGITAL 

STAGE 

FILTER COEFFICIENTS 
VALUE 

1 

0. 10000000D 

02 

2 

-0 • 1.6666665D 

02 

3 

0. 133 333 33 D 

02 

4 

-0 . 10000000D 

02 

5 

0 .66666671 D 

01 

6 

-0.33333345D 

01 


Fig. C-2. LINPO Output for Raised-Cosine Filter Design. 
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C.1.3 Flowchart and Program Listing 

A flowchart and program listing for the C0CAL1 subroutine are 
given in Figs. C-3 and C-4. The left-hand half of the flowchart 
consists of two loops which calculate coefficients required for the 
LINPO A array. The outer loop (thru statement 210) runs from 1 to 
NEQ with the first execution corresponding to the objective function 
and the remaining NEQ-1 executions corresponding to the constraint 
equations. For each execution of the outer loop, there are NVA 
executions of the inner loop (thru statement 208) to calculate the 
NVA coefficients of each equation. That set of NVA coefficients 
forms one row of the A array. Since the objective function will 
normally be a combination of equations, the appropriate coefficients 
of the equations that make up the objective function must be combined 
during the first passage through the outer loop. This combining of 
coefficients is controlled by the NOBJ array. 

The right-hand half of the flowchart represents three operations. 
First, the coefficients in the first row of the A array are multiplied 
by -1 if the objective function is of the type (C.3c) or (C.3d). As 
stated in Table B.l, the objective-function coefficients must be multi- 
plied by -1 before a linear-programming solution is obtained, assuming 
that the objective-function value Z is to be maximized, which always 
occurs internally in the linear-programming subroutine. However, if 
the coefficient signs are not changed the value -Z is maximized and is 
actually a minimization of Z. Consequently, if Z is to be minimized 
the signs are not changed and by default a minimization problem is 
changed to a maximization problem. Next, the constraint equations 
must be checked to see if the right-hand sides are negative. If they 



are, both sides of the equation are multiplied by -1 and the ine- 
quality changed. This results in all non-negative numbers in its 
B array. Finally, LINPO is called to solve the problem which has 
now been properly formulated. 




Fig. 03. Flowchart for Subroutine C0CAL1. 
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C 

C 

C 

C 


C 


200 

201 

202 


C 

c 


203 


204 

205 

206 


208 

209 

210 
C 

C 


SUBROUTINE C0CAL1 (T,NTYPE,B2 , IE ,NOBJ,M, T 1 , NEO , NV A , N SP , NT ) 

THIS SUBROUTINE CALCULATES AND SETS UP ARRAYS REOUIRFD 
FOR THE L IN PO DESIGN OF FREQUENCY-SAMPLING FILTERS. 

DIMENSION NTYPE(NEQ)»IE(NEQ)» NOB J ( 1 ) 

REAL* 8 SUM,B1,B2(NE0) ,DCOS , P I , T ( 1 ) , A ( 27 , 50 ) , M , T1 , 

IK, TT,DSIN 
P 1 = 3. 141 5926D0 

CALCULATE COEFFICIENTS FOR ARRAY A (THRU STATEMENT 210) 

DO 210 1=1, NEO 
NT YP=NTYP E ( 1 ) + l 
GO TO ( 200,201 ,200,201 ) ,NTYP 
NT = 0 

GO TO 202 
NT = 1 
K =0 .0 DO 
IS GN = 1 
I 1 1 = 1 

IF (I.EO.l) I 1 1 = NO B J ( 1 ) 

CALCULATE COEFFICIENTS FOR ROW I OF ARRAY A 
(Thru statement 208) 

DO 208 J = 1 , NVA 
SUM=0. ODO 
DO 206 I 1 = 1,11 I 
JJ=I-1 

I F ( I .GT.l ) GO TO 203 
JJ=II+1 

N BS= I AB S ( NOB J ( JJ ) ) 

ISGN = ISIGN( 1,N0BJ( JJ) ) 

J J=NBS 

TT = T ( J J ) / T1 

B1=(2.0D0*PI*K )/M 

IF ( NT -0)204,204,205 

SUM=SUM+DC0S(B1*TT)*ISGN 

GO TO 206 

SUM=SUM-B1*DS I N ( B 1*TT )*I SGN 
CONTINUE 

A ( I , J )=SUM*(2.0D0/M) 

IF (K.EO.O.ODO) A { I ,J)=A(I ,J)/2.0D0 

K =K + 1 ,0 DO 

CONTINUE 

DO 209 J =2 , NVA , 2 
A ( I , J ) =-A ( I , J ) 

CONTINUE 

CONVERT MINIMIZATION PROBLEM TO A MAXIMIZATION PROBLEM 
(THRU STATEMENT 212) 

NTYP=NTYPE( 1 )+l 


Fig. C-4. Program Listing for Subroutine C0CAL1. 
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GO TO (211 ,211 ,213, 213), NTYP 

211 DO 212 J = 1 , NVA 
A ( 1 , J) =-A ( 1 , J ) 

212 CONTINUE 

C CONVERT CONSTRAINT EQUATIONS IF RIGHT-HANDED SIDES ARE 

C NEGATIVE (THRU STATEMENT 218) 

213 DO 218 1=2, NEO 

I F( 82 ( I ) 1214,218,218 

214 R2 ( I ) = -B2 ( I ) 

DO 215 J = 1 , NVA 

215 A ( I , J ) =-A ( I , J ) 

I F ( I E ( I )— 2 1216,218,217 

216 I E ( I 1 = 3 
GO TO 218 

217 I E ( I 1 = 1 

218 CONTINUE 

CALL LINPO(A,B2, IE , NEQ , NV A , NS P ,NT 1 

RETURN 

END 


Fig. C-4 (Continued). Program Listing for Subroutine COCAL1. 
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C. 2 SUBROUTINE COCAL2 


COCAL2 is a double-precision computer subroutine that formulates 
an equalizer design problem for a linear-programming solution in the 
same manner that C0CAL1 structures a frequency-sampling filter design 
problem. The expression which was derived in Section 4.1 for the out- 
put of a data channel equalized with a frequency-sampling filter was 

H L 

y(t) = ^ V c(t - mT) + 
m=0 


Z l 

m =n v= i 


(-D 


, 2H. 
k k 

N 


cos (2irkm/N) c(t — mT) 


(C. 4) 


and similarly for a transversal equalizer 

L 

y (t) = Y' (-1)“ H c(t - mT) . (C.5) 

/ . m 

m=0 


It was also shown that the c(t) values must be obtained from the actual 
channel impulse response by 


c(t) 



c (x) dx 
c 


c (r) dx 
c 


0 < t < T 


T < t < ® , 


(C. 6) 


and that the derivatives of (C.4) and (C.5) are formed by replacing 
the appropriate c'(t - mT) terms with 
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c (t) 
c 


0 < t < T 


c'(t) =< 


c (t) - c (t - T) t > T 
c c 


(C. 7) 


During a LINPO design of either type of equalizer, the coefficients of 

the H or H terms must be calculated from (C.4) or (C.5) along with 
k m 

(C. 6) and (C.7). In addition, other LINPO input arrays must be gener- 
ated. 


C.2.1 User Instructions 

The call statement for the subroutine is CALL C0CAL2 (CT,NTYPE,B, 

IE ,NOBJ ,N ,T ,NEQ ,NVA,NSP ,NT,CC) where, with the following exceptions, 
the input parameters have been defined previously in Tables B.l and 
C.l. Acceptable forms of the objective function are still defined 
by (C.3), with h(t^) replaced by y(t__). The possible values in the 
array NTYPE must be expanded to include both frequency-sampling' and 
transversal filters. Values of NTYPE from 0 to 3, which previously 
identified the type of objective function and constraint equations 
for frequency-sampling pulse-shaping filters, now apply to frequency- 
sampling digital filter equalizers with the h(t) in Table C.l replaced 
by y (t) . NTYPE values from 4 through 7 identify transversal-equalizer 
designs and are correspondingly defined if each integer in the NTYPE 
description is increased by 4, e.g., 0 by 4, 1 by 5, etc. As an 
example, if the objective function for a transversal equalizer design 
is in the form of (C.3c), i.e., the minimization of a combination of 
y(t^) equations, the first value in the array would be NTYPE = 6. 

Unlike frequency-sampling filter designs, the number of unknown 
multiplier coefficients for a transversal design is equal to the number 
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of impulse- response values. Hence, NVA must be set equal to N for 
LINPO designs of transversal filters. 

The parameter CC is not an input value but represents an external 
function that is called by C0CAL2. Consequently, a function subprogram 
CC, which will be described in Section C.4, must be used with C0CAL2. 

C.2.2 Example Problem 

An equalizer design problem introduced in Section 4.2.1 was to 

determine both frequency-sampling and transversal filters to equalize 

2 2 7T t 

the channel c (t) = 4ir t e u(t) to meet the set of constraints 

c 

specified by (4.16) and the objective function of (4.17). Filters 
with a sampling interval of 0.1 seconds and 25 impulse response 
values were specified. For the frequency-sampling filter, 12 multi- 
plier coefficients were specified and from the impulse-response speci- 
fication, 25 coefficients were allowable for the transversal filter. 

The input program required for a C0CAL2 solution of this problem is 
illustrated in Fig. C-5. 

Only one main program was required for this example since the set 
of specifications on the desired pulse shape were the same for both 
equalizers. C0CAL2 is called twice, first for a frequency-sampling 
design and next for a transversal design. The only input-data dif- 
ference between the two runs are the set of values for NTYPE and the 
value of NVA. Subroutine CHECK which is also called for each run and 
its associated input parameters will be described in Section C.3. 

The statement NN=2 specifies which channel model in subprogram CC is 
being equalized. LINPO outputs for the two runs are illustrated in 
Figs. C-6 and C-7. As indicated, the frequency-sampling filter 
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requires 11 multipliers while the transversal filter only requires 
three. 


DIMENSION NTYPE(13) , IE(13) ,NOBJ (5) 

REAL*8 CT(12) ,B(13) ,N,T,TIN,TM,DT 
COMMON NN 
EXTERNAL CC 

DATA CT/0. 5DO,0 . 75DO, 1. DO, 1. 5DO, 2 . DO,2 , 25DO ,2. 5D0, 3. DO 
1,0. 5DO, 1. DO ,2 . DO , 2 . 5D0/ 

DATA NTYPE/1, 8*0, 4*1/ 

DATA B/4*O.DO,1.DO,8*O.DO/ 

DATA IE/0 ,8*2 ,2*3, 2*1/ 

DATA NOBJ/4,-1,-3,5,7/ 

DATA N,T ,NEQ,NVA,NSP/25 . D0,0. IDO, 13, 12,0/ 

NN=2 

DATA TIN ,TM,DT/0 . DO , 3. DO ,0 . 05D0/ 

CALL COCAL2 (CT ,NTYPE , B , IE ,NOBJ ,N ,T ,NEQ ,NVA,NSP ,NT , CC) 
CALL CHECK(TIN,TM,DT,N,T,NVA,NT,CC) 

DO 1 1=1,13 
1 NTYPE ( I ) =NTYPE ( I ) +4 
NVA=25 

CALL COCAL2 (CT , NTYPE , B , IE ,NOB J ,N , T ,NE0 ,NVA,NSP ,NT , CC) 
CALL CHECK (TIN , TM , DT , N , T , NVA , NT , CC) 

STOP 

END 


Fig. C-5. Main Program for COCAL2 Solution. 

C.2.3 Flowchart and Program Listing 

Since COCAL2 performs the same functions using y(t) specifica- 
tions as COCAL1 does with h(t) specifications, the flowchart of Fig. 
C-3 is still valid for C0CAL2, with the exception that the flowchart 
statement numbers 208, 210, and 218 are now 316, 322, and 330, re- 
spectively. Also, coefficients in y(t^) and y'(t^) are calculated 
instead of h'(t^) and h(t^). The COCAL2 program listing is furnished 
in Fig. C-8. 
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OAT A LOADFD 
FEASIBLE 


OPTIiiUp Z 0.0 


OASIC SOLUTION 


VAR I A B L F 

VALUE 


: 1 

0.175575220 

0 1 

2 

0.219808290 

01 

3 

0.387010080 

0 1 

5 

0.891418840 

01 

4 

0.423829570 

01 

7 

0.934511220 

01 

R 

0.887938430 

0 1 

o 

0.820470710 

01 

10 

0.724330660 

01 

A 

0. 982218 71/) 

01 

1 1 

0.366461890 

01 

1 6 

0 . 0 



DIGITAL 

FILTER COEFFICIEl 

STAGE 

VALUE 


1 

0.175575220 

01 

2 

-0.439616570 

01 

■ 3 

0.774020160 

01 

4 

-0.124765910 

02 

5 

0.178283770 

02 

6 

-O’. 1.96443740 

02 

7 

0. 186902240 

02 

8 

-0.177587690 

0 2 

o 

0. 164094140 

02 

10 

-0.144866130 

0 2 

1 1 

0.732923790 

01 

12 

0.0 



Fi 8* C 6 * LINPO Output for Frequency-Sampling Equalizer Design. 



DATA LOADFD 
FEASIBLF 

OPTIMUM Z 0.0 


BASIC SOLUTION 


VARIABLE 

VALUE 


5 

0.0 


4 

0.0 


9 

0.0 


23 

0.0 


R 

0.0 


13 

0.48966586D 

01 

19 

0.5 64 4405 8 D 

00 

0 

0.0 


18 

0. 12696229D 

01 

7 

0.0 


22 

0.0 


25 

0.0 



DIGITAL 

FILTER COEFFICIENTS 

STAGE 

VALUE 

1 

0.0 

2 

0.0 

3 

0.0 

4 

0.0 

5 

0.0 

6 

0.0 

7 

0.0 

8 

0.0 

9 

0.0 

10 

0.0 

11 

0.0 

12 

0.0 

13 

0.48966586D 01 

14 

0.0 ; ' 

15 

0.0 

16 

0.0 

17 

0.0 

18 

-0.126962290 01 

19 

0.564440580 00 

20 

0.0 . 

21 

0.0 

22 

0.0 

23 

0.0 

24 

0.0 

25 

0.0 


Fig. C-7. LINPO Output for Transversal Equalizer Design. 
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SUBROUTINE COCAL2 <T,NTYPE,B2, IE,NOBJ,M, T1 , NEO , NVA , M SP , NT 
1 ,CC ) 

C 

C THIS SUBROUTINE CALCULATES AND SETS UP ARRAYS REOUIR ED 

C FOR THE LINPO DESIGN OF FREQUENCY-SAMPLING OR 

C TRANSVERSAL EQUALIZERS. 

C 

DI MENS ION NTYPE (NEO ) , I E ( N EO ) , NOB J { 1 ) 

REAL*8 SUM » B2 ( NEO ) ,CC » DCOS , P I » D , T { 1 ) , A ( 27 , 50 ) » Y» M » 
1T1,K,N,X,Z 
R EAL* 8 A 1 { 2 5 , AO ) 

EXTERNAL CC 
P 1=3. 1A1 5926D0 

C CALCULATE COEFFICIENTS FOR ARRAY A (THRU STATEMENT 322) 

DO 322 1=1, NEO 
NT YP=N TYP E ( I ) + l 

GO TO (300,301,300,301,302,303,302,303) ,NTYP 

300 NT =0 

GO TO 30 A 

301 NT = 1 

GO TO 30 A 

302 NT = A 

GO TO 30 A 

303 NT =5 

30A K =0 .0 DO 

ISGN = 1 
I 11 = 1 

IF (I.EO.l) III=N0BJ(1 ) 

N VA1=NVA 

IF(NT.GE.A) NVA 1=1 

C CALCULATE COEFFICIENTS OF ROW I OF ARRAY A 

C (THRU STATEMENT 316) 

DO 316 J = 1 , NVA 1 
SUM=0. ODO 
DO 313 I 1=1, I I I 
JJ = I-1 

I F ( I .GT.l ) GO TO 305 
JJ = I 1 + 1 

N BS= I AB S ( NOB J ( JJ ) ) 

ISGN = ISIGN( 1, NOB J ( J J ) ) 

J J=NB S 

305 L= T ( J J ) /T 1 + 1 

TTL=FLOAT(L) 

IF (DBLE(TTL).GT.M) L=M 
N=0 .ODO 
DO 311 I J =1 , L 
X=T(JJ)-N*T1 
Z= X-Tl 


Fig. C-8. Program Listing for Subroutine C0CAL2. 
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IF U) 306,306,307 

306 I F(NT.E0.0.0R.NT.EQ.4) CALL DOG 16 ( O.ODO , X , CC , Y ) 

IF (NT. EO. 1 .0R.NT.EQ.5 ) Y=CC(X) 

GO TO 308 

307 IF (NT.E0.0.0R.NT.EQ.4) CALL DQG1 6 < Z , X , CC , Y ) 

I F( NT.EO.l .OR. NT. EO. 3 ) Y=CC ( X ) -CC ( Z ) 

308 IF ( NT . GE . 4 ) GO TO 309 
D=(2.0D0*PI*K*N)/M 
SUM=ISGN*Y*DCOS ( D) +SUM 
GO TO 310 

309 A1 ( I I, IJ)=Y*ISGN 

310 N=N+1.0D0 

311 CONTINUE 

I F(NT.LT.4.0R.L.E0.NVA) GO TO 313 
LL = L + 1 

DO 312 J I=LL , NVA 

312 A1 ( I I , J I ) =0 . ODO 

313 CONTINUE 

I F ( NT. GE . 4 ) GO TO 316 
A ( I , J )=SUM*(2.0D0/M) 

314 IF(K.EO. O.ODO) A ( I , J ) = A ( I , J ) / 2 .ODO 

315 K =K + 1 .ODO 

316 CONTINUE 

I F (NT.LT.4) GO TO 320 

317 DO 319 I J =1 »NVA 
S UM=0 .ODO 

DO 318 11=1, III 

318 SUM=SUM+A1 ( I I, I J ) 

319 A( I , I J ) = S UM 

320 CONTINUE 

DO 321 J=2 , NVA , 2 

321 A ( I , J )=-A( I, J ) 

322 CONTINUE 

C CONVERT MINIMIZATION PROBLEM TO A MAXIMIZATION PROBLEM 

C (THRU STATEMENT 324) 

N TY P= NT Y PE ( 1 ) + 1 

GO T0( 323,323,325,325,323,323,325,325) ,NTYP 

323 DO 324 J=1,NVA 

A ( 1 , J ) =— A ( 1 , J ) 

324 CONTINUE 

C CONVERT CONSTRAINT EQUATIONS IF RIGHT-HANDED SIDES ARE 

C NEGATIVE (THRU STATEMENT 330) 

325 DO 330 1=2, NEO 

I F ( B 2 ( I ) )326,330,330 

326 B2 ( I ) = — B2 ( I ) 

DO 327 J = 1 , N VA 

327 A ( I , J ) =-A ( I , J ) 

I F( IE ( I )-2 ) 32 8 ,330,329 


Fig. C-8 (Continued). Program Listing for Subroutine C0CAL2. 
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328 T EC I ) =3 
GO TO 330 

329 I E { I ) =1 

330 CONTINUE 

CALL LINP0(A,B2, I E ,NEQ,NVA , NSP , NT ) 

RETURN 

END 


Fig. C-8 (Continued). Program Listing for Subroutine COCAL2 . 
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C. 3 SUBROUTINE CHECK 

Equations (C.4) and (C.5) represent the received pulse at the 
output of a channel equalized with a frequency-sampling filter and 
transversal filter respectively. In order to verify the results of 
the LINPO equalizer designs, these equations along with (C.6) were 
programmed to form the double-precision computer subroutine called 
CHECK which calculates and generates plots of y(t). 

The call statement for the subroutine is CALL CHECK (TIN, TM,DT, 

N ,T ,NVA,NT, CC) where the input parameters TIN, TM, and DT are defined 
in Table C.2 and the other parameters have been defined previously. 

TABLE C.2 

CHECK INPUT PARAMETERS 

TIN: A double-precision constant equal to the initial time 

point at which y(t) is desired to be evaluated. 

TM: A double-precision constant equal to the maximum time 

point at which y(t) is desired to be evaluated. 

DT: A double-precision constant equal to the time interval 

desired between consecutive evaluations of y(t). 

Filter coefficients or H^ need not be externally entered into the 
subroutine since these values are placed in a labeled COMMON block by 
LINPO at the completion of the design and are therefore available to 
CHECK. The integer constant NT is required by CHECK to determine which 
equalizer has been designed, but is generated by C0CAL2 and does not 
have to be stated in the main program. An example of the use of the 
CHECK subroutine in a main program is illustrated in Fig. C-5 with 
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the output plots generated by the subroutine given in Figs. 4-3 (b) 
and 4-4 (b) . 

A flowchart and program listing for the subroutine are presented 
in Figs. C-9 and C-10 respectively. The primary loop of the program 
is executed for each of the NDATA time points at which the function 
y(t) is to be evaluated. During each execution, channel impulse- 
response values obtained by integrating c^Ct) from the function sub- 
program CC are combined with either the or values from LINPO, 
depending on the values of NT, to form y(t^). After all NDATA values 
of y(t^) have been calculated, GRAPH1 [52] is called to plot y(t^) 
versus t . 
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Fig. C-9 . Flowchart for Subroutine CHECK. 
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SUBROUTINE CHECK < T IN , TM , DT , M , T 1 ,NVA,NT,CC) 

C 

C THIS SUBROUTINE CALCULATES THE IMPULSE RESPONSE OF 

C AN EQUALIZED DATA TRANSMISSION CHANNEL. EITHER 

C TRANSVERSAL OR FREQUENCY-SAMPLING DIGITAL FILTERS 

C MAY BE US EO AS EQUALIZERS. 

C 

DOUBLE PRECISION H ( 40 ) 

DIMENSION A( 500) , YY( 500 ) 

REAL*8 SUM,CC,DC0S,PI,D,V,Y(200) , M , T 1 , K , N , X , Z ,SUM 1 

REAL*8 TIN,TM,DT,SUM2,TMP 

COMMON/HNAME/H 

EXTERNAL CC 

PI =3 . 1 41 592 6D0 

T MP = T IN 

NDATA= (TM-T IN ) /DT+1 

C CALCULATE VALUE OF Y ( T ) FOR EACH TIME T (THRU 

C STATEMENT 110) 

DO 110 I=1,NDATA 

SUM2=0.0D0 

A ( I )=SNGL( TIN) 

L=TIN/T1 + 1 
TTL=FLOAT( L) 

IF (DBLE(TTL).GT.M) L=M 
N =0.0 DO 

C CALCULATE IMPULSE RESPONSE VALUES OF UNEOUALIZED 

C CHANNEL (THRU STATEMENT 104) 

DO 104 IJ=1,L 
X=TIN-N*T1 
Z= X-Tl 

IF (Z) 100,100,101 

100 CALL DQG16(0.0D0,X,CC,V) 

GO TO 102 

101 CALL DQG16( Z,X,CC,V) 

102 N=N+1.0DO 
IF(NT.LT.4) GO TO 103 

I F( H( IJ ) .EQ.O.ODO ) GO TO 104 
SUM2=SUM2+V*H( I J ) 

GO TO 104 

103 Y( I J ) = V 

104 CONTINUE 

C DETERMINE Y ( T ) FOR TRANSVERSAL FILTER 

I F( NT.GE.4 ) YY( I ) = SUM2 
IF (NT.GE.4) GO TO 109 
S UM1 = 0. ODO 

105 K= 0 . ODO 

C CALCULATE Y ( T ) FOR FREQUENCY-SAMPLING FILTER 

C (THRU STATEMENT 109) 


Fig. C-10. Program Listing for Subroutine CHECK. 



no 108 J = 1 * IM VA 

IF (H(J).EQ.O.ODO) GO TO 107 

SUM=0 .000 

N= 0 • 0 D 0 

DO 106 I J= 1 » L 

0= ( 2 .0 DO* P I #K*N ) /M 

S UM= Y ( I J ) *DCOS ( 0 ) +SUM 

N= N+ 1 . ODO 

106 CONTINUE 

SUM1 =S UM1 +SUM*H ( J ) 

107 K =K + 1 .000 

108 CONTINUE 

Y Y ( I ) =SNGL (SUM1/M ) 

109 TI N=T I N+DT 

110 CONTINUE 
TI N=TMP 

CALL GRAPH1(A,YY,NDATA,1 ) 
.RETURN 
END 


Fig. C-10 (Continued). Program Listing for Subroutine CHECK. 
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0.4 FUNCTION SUBPROGRAM CC 


Function subprogram CC is called by the subroutines COCAL2 and 
CHECK to calculate channel impulse-response values. The eight models 
of typical channels which are available from the subprogram are listed 
in Table C.3. It should be noted that the channels have been normal- 
ized in the subprogram to have a unity breakpoint frequency. 


NN 

1 

2 

3 

4 


TABLE C. 3 

IMPULSE-RESPONSE CHANNEL MODELS 


c (t) 
c 


0 -2irt ,. N 
2ir e u(t) 


NN , C c (t) 

5 2^(2ir) 5 t 4 e 2lTt u(t) 


2 -2irt * 
it te u(t) 


1 .6 5 -2irt ... 

(2tt) t e u ( t ) 


120 


1 ,„ \ 3 2 -2irt . s 

2 -( 2 tt) t e u(t) 


7 e 2lTt sin (2irt) u(t) 


l/o _3 -2irt . . 
■g(2Tr) t e u(t) 


8 e 2TTt cos ( 2 tt t ) u(t) 


In order to specify which channel model is desired, the parameter 
NN must be declared with a blank COMMON statement and its value given 
in the main program that calls either C0CAL2 or CHECK. The main program 
must also declare the function CC to be external. For example, see 
Fig. C-5 where the statement NN = 2 specifies the second channel re- 
sponse listed in Table C.3. The CC function program listing is given 
in Fig. C-ll. 
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DOUBLE PRECISION FUNCTION CC(T) 

C 

C THIS FUNCTION SUBPROGRAM EVALUATES ONE OF EIGHT 

C POSSIBLE IMPULSE RESPONSE FUNCTIONS FOR THE 

C TIME POINT T. 

C 

DOUBLE PRECISION T ,DE XP , AL PHA ,ALPN , AT , TN ,PROD, DS I N , DCOS 

COMMON NN 

N=NN 

A LPHA =6 .283185200 
AT=ALPHA*T 

IF (AT. GE. 165. DO) GO TO 601 
IF(N.LE.6) GO TO 603 
M=N— 6 

GO T0< 604,605) ,M 

603 I F ( N. EQ . 1 ) GO TO 607 

IF (T.LT.l .D-9) GO TO 601 
TN=T**(N-1) 

ALPN=ALPHA**N 
GO TO 608 

607 TN = 1.0DO 

A LPN= 1 . ODO 

608 PR0D = 1.0D0 
DO 600 1=2, N 

600 PROD=( 1-1 )*PROD 

C CALCULATE CC ( T ) FOR CHANNEL WITH SIMPLE POLES 

CC= ( DE XP ( -AT ) *AL PN*TN ) /PROD 
GO TO 602 

C CALCULATE CC(T) FOR CHANNEL WITH COMPLEX POLES 

604 CC=DEXP(-AT)*DSIN(AT) 

GO TO 602 

C CALCULATE CC ( T ) FOR CHANNEL WITH COMPLEX POLES 

C AND A SIMPLE ZERO 

605 C C=DE XP ( — A T ) *DCOS ( AT ) 

GO TO 602 

601 C C=0 • ODO 
GO TO 602 

606 CC=1 • ODO 

602 . RETURN 

END 


Fig. C-ll. Program Listing for Function CC. 
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